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FOREWORD 


The Computational Fluid Dynamics (CFD) computer codes and Knowledge-Based System 
(KBS) were generated under NASA contract NAS3-25644 originating from the Office of 
Advanced Concepts and Technology and administered through NASA-Lewis Research Center. 
The support of the Program Manager, Anita Liang, and the advice and direction of the 
Technical Monitor, Robert Hendricks, are gratefully appreciated. Major contributors to code 
development were: 

• Dr. Bharat Aggarwal: KBS and OS/2 PC conversion of labyrinth seal code KTK 

• Dr. Antonio Artiles: cylindrical and face seal codes ICYL and IFACE 

• Dr. Mahesh Athavale and Dr. Andrzej Przekwas: CFD code SCISEAL 

• Mr. Wilbur Shapiro: gas cylindrical and face seal codes GCYLT, GFACE, and seal 
dynamics code DYSEAL 

• Dr. Jed Walowit: spiral groove gas and liquid cylindrical and face seal codes SPIRALG 
and SPIRALI. 

The labyrinth seal code, KTK, was developed by Allison Gas Turbine Division of General 
Motors Corporation for the Aero Propulsion Laboratory, Air Force Wright Aeronautical 
Laboratories, Wright-Patterson Air Force Base, Ohio. It is included as part of the CFD 
industrial codes package by the permission of the Air Force. 
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NOMENCLATURE FOR CODE SPIRALG 


A 

A, 

{a} 

o k 

a Lj 

B xy 


BB xy 

B<t>v 


BB, X 

B 0 

[B] 

{b} 

C 

[C j ] 

[C J ] 

[C j ’ k ] 

[O] 

[C* j ] 

[V] 

[D*] 

[I> k ] 


= dimensionless flow control area, (area/R 0 2 ) 

= coefficients of second-order linear operator defined by Eq. (2-44), 1=1,.. .,6 
= column vector of eccentricity coefficients defined by Eq. (2-41) 

= kth component of {a} evaluated at grid point (i,j) 

= damping coefficient relating force in x direction to velocity in y direction, B xx , 
B yx , B yy and B H are similarly defined 

= damping coefficient relating moment about x axis to angular velocity about y 
axis, B w , B^ and B w are similarly defined 

= damping coefficient relating force in x direction to angular velocity about x 
axis, B xv , B yV)/ , B^ and B ZV) , are similarly defined 

= damping coefficient relating moment about x axis to velocity in x direction, 

B,,*, B^ y , B w , B^ and B^ are similarly defined 

= dimensionless damping coefficient B xy /B 0 , same definitions apply to B xx , B^, B yy , B 

= dimensionless damping coefficient B^BqRq 2 ), same definitions apply to B^, B^, 
B w 

= dimensionless damping coefficient B^BqRq), same definitions apply to B xy , 

B y( j>, B yX p B Z( j,, B z¥ 

= dimensionless damping coefficient B^AB^), same definitions apply to B 

B<j>y’ B\jry? B^, B^ 

= characteristic damping constant, 12 jliR 0 4 /C 3 
= dimensionless damping coefficients in matrix form 

= column vector of coefficients of e', arising from linearization of Eq. (2-39) 

= clearance (cylindrical seal) or reference film thickness (face seal) 

= coefficient matrix used in Newton-Raphson linearization procedure, see Eq. (2-34) 
= coefficient matrix obtained from steady state solution 
= derivative of [C J ] with respect to £ k 

= diagonal coefficient matrix whose components are given by Eq. (2-51) 

= complex coefficient matrix used in complex stiffness solution, [C j ] + 3 g[C j ] 

= coefficient matrix used in Newton-Raphson linearization procedure, see Eq. (2-34) 
= coefficient matrix obtained from steady state solution 
= derivative of [!>] with respect to e k 
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[E] 

[E j ] 

[E k ] 


e x’ e y 


Fy 

G 

H 

H r 

h g 

h r 

ij»k 

y 

3 

K xy 


K 


w 


K 


X<|) 


k* k 


K 

K 

K 

Ko 

[KJ 

K 

L 


= coefficient matrix used in Newton-Raphson linearization procedure, see Eq. (2-34) 
= coefficient matrix obtained from steady state solution 
= derivative of [E J ] with respect to e k 
= eccentricity in x,y direction (cylindrical seal) 

= axial displacement, (face seal) 

= residual outflow function from flow balance at grid point (i,j) 

= second-order nonlinear operator defined by Eq. (2-39) 

= steady state, unperturbed, value of H r 
= dimensionless film thickness, h/C 
= film thickness over grooves, see Fig. 2 
= film thickness over ridges, see Fig. 2 
= subscripts used generically as indices 
= unit vectors in 0,s directions 
= unit imaginary number, V-l 

= stiffness coefficient relating force in x direction to displacement in y direction, 
K xx , K yx , K yy and are similarly defined 

= stiffness coefficient relating moment about x axis to rotation about y axis, K^, 
and K w are similarly defined 

= stiffness coefficient relating force in x direction to rotation about x axis, K xr 
K y(j) , K yxp and are similarly defined 

= stiffness coefficient relating moment about x axis to displacement in x 
direction, K^, K^, K w , K^ z and are similarly defined 

= dimensionless stiffness coefficient K^/K^, same definitions apply to K xx , K^, K^, 

= dimensionless stiffness coefficient K^K^ 2 ), same definitions apply to K^, K^, 

= dimensionless stiffness coefficient K X( /(K 0 R 0 ), same definitions apply to K xvu , 

i> t> K K ^ 

- rs 'y<)>’ 

= dimensionless stiffness coefficient ^/(KqRq), same definitions apply to K^, 

Rv|/y» K<j)zi K\|/ Z 

= characteristic stiffness constant, p 0 Ro 2 /C 
= dimensionless stiffnesses in matrix form 
= spiral groove coefficient defined by Eq. (2-25), i=l,2,...,8 
= seal length, see Fig. 1 
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x 



L 

l g 

It 

££ 

M 

M x ,M y 

M x ,M y 

N 

N . 

rt 

n p 

P 

P 

Py 

Pi 

P,Pr 

{P'l 

{p; k i 

{p*} 

{p* k } 

{p*i 

{p 3 } 

p 

p' 

Po 


Pl’Pr 

0 

Qs’Qe 

Qtj’Qii 

Qin 


= dimensionless length, L/(2R 0 ) 

= groove width, rA0 g 
= ridge width, rA0 r 

= second-order linear operator defined by Eq. (2-44) 

= number of grid points in s direction 
= applied moment about x,y axis 
= dimensionless moment, (M x ,M y )/(R 0 3 p 0 ) 

= number of grid points in 0 direction 
= number of spiral grooves 
= unit vector normal to S 
= unit vector normal to groove 
= dimensionless pressure, (p-p 0 )/p 0 
= steady state unperturbed value of P 
= dimensionless pressure, P, at grid point (i,j) 

= dimensionless pressure, P, at point i shown in Fig. 3, i=l,...,9 
= dimensionless boundary pressures (prP 0 )/po ? (PrPoVPo 

= column vector of dimensionless pressure disturbances due to perturbation in {e} 
= column vector of disturbance pressures associated with e k 
= complex amplitude of {P'J, {P # } = {P*}e 3c ” 

= column vector of complex stiffness pressures associated with e k 
= real part of {P*} 

= imaginary part of {P*} 

= global pressure in absolute units 
= local pressure in absolute units 

= reference pressure in absolute units, normally taken to be the minimum of the 
boundary pressures 

= boundary pressures in absolute units at s,,s r 
= dimensionless flow vector, 12pR 0 ^/(p 0 C 3 ) 

= components of dimensionless flow vector in s,0 directions 
= dimensionless flow components shown in Fig. 3 
= dimensionless flow rate, 1 2pq in /(p 0 C 3 ) 
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xi 





q s >q 9 

qA 

q in 

4' 


q^qe 

qX 

R 

Ro 

Ri 

{R J } 

{R J } 

{R j,k } 

{R{’ k } 

{R j,k } 

r 

S 

Si,S r 

s 


Si 

s r 

T 

T 

T c 

T f 

T r 


- global flow vector, mass flow rate per unit transverse length divided by density 
at pressure p 0 

= components of q in s,0 directions 

= global mass flow rate per unit area displaced by rate of decrease of film, 
divided by density at p 0 

= volumetric flow rate measured at pressure p 0 

= local flow vector, mass flow rate per unit transverse length divided by density 
at pressure p 0 

= components of q' in s,0 directions 

= local mass flow rate per unit area displaced by rate of decrease of film, divided 
by density at p 0 

= dimensionless coordinate, r/R 0 , taken as 1 for cylindrical seal 

= reference radius, taken as outside radius for face seal and shaft radius for 
cylindrical seal 

= inside radius of face seal 

= column vector used in Newton-Raphson linearization procedure, see Eq. (2-34) 
= column vector obtained from steady state solution 
= derivative of {R j } with respect to e k 
= column vectors whose components are given by Eq. (2-52) 

= complex column vectors used complex stiffness solution {R j,k } - 3cj{R j,k } 

= radial coordinate, taken as R 0 for cylindrical seal 
= dimensionless coordinate, s/R 0 
= dimensionless boundary coordinates s/R 0 , s/R 0 
= dimensionless length of line surrounding flow control area (length/R 0 ) 

= transverse coordinate, s = r for a face seal and s = z for a cylindrical seal 
= transverse coordinate at start of groove 

= left boundary coordinate, Sj = R ; for face seal and s t = -L/2 for cylindrical seal 
= right boundary coordinate, s r = R 0 for face seal and s r = L/2 for cylindrical seal 
= torque 

= dimensionless torque, T/(p 0 R 0 2 C) 

= T obtained with course grid in Romberg extrapolation example 
= T obtained with fine grid in Romberg extrapolation example 
= T obtained by Romberg extrapolation 


xii 
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t 

t 

b 

Ul 

u 2 

W x ,w y 

W x ,w y 

W z 

W z 

{W} 

{W} old 

x,y 

z 


^opt 

p 

Pop, 

p 

r 

aA 

aA s 

Ap 

a p ; 

ap; 

AStj 

ASjj 

As 

A0 


= time 

= dimensionless time, cot/(2A) 

= unit vector tangential to groove 
= velocity of grooved surface 
= velocity of smooth surface 
= applied loads in x,y direction (cylindrical seal) 

= dimensionless loads (W x ,W y )/(p 0 R 0 2 ) 

= applied load in z direction (face seal) 

= dimensionless load, W^p^ 2 ) 

= column vector containing dimensionless loads and moments {W x ,W y ,M x ,M } 
for cylindrical seal, {W z ,M x ,M y } for face seal 

= {W} at previous iteration in eccentricity homing process 

= coordinate variables, see Fig. 1 

= dimensionless axial coordinate for cylindrical seal, z/Rq 

= axial coordinate, measured from axial center for cyl. seal or from reference 
film, C, for face seal 

= groove to pitch ratio, l g /(l g +l r ) 

= value of a for maximum stagnation pressure gradient in concentric cylindrical seal 
= spiral groove angle, angle between grooves and surface velocity 
- value of p for maximum stagnation pressure gradient in concentric cylindrical seal 
= numerical damping factor used in eccentricity homing process 
= film thickness ratio, hg/h,. 

= dimensionless flow control area about single grid point, shaded area in Fig. 3 
= portion of AA in quadrant containing point i, see Fig. 3 
= global pressure difference, see Fig. 2 
= pressure difference across groove, see Fig. 2 
= pressure difference across ridge, see Fig. 2 
= line or arc length associated with Q}j 
= line or arc length associated with Qjj 
= transverse length of groove-ridge pair 

= circumferential extent of groove-ridge pair, see Fig. 2; (also used generally for 
change in 0) 
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A e s 

Ae r 

8 


£ £ 
0 X ,O y 


e z 

e k 

[e] 

{£} 

{ £ }old 

£' 

n 

e 

e 

A 

A§ 

F 

a 

x 

X 

x' 

O 

¥ 

Q 

Q 

co 

CO 

«! 


= circumferential extent of groove, see Fig. 2 
= circumferential extent of ridge, see Fig. 2 
= dimensionless groove depth, (h -h r )/C 

= value of $ for maximum stagnation pressure gradient in concentric cylindrical seal 
= eccentricity ratio, (e x ,e y )/C (cylindrical seal) 

= axial displacement ratio eJC, (face seal) 

= kth component of eccentricity matrix 

= row matrix of eccentricity components, [£ z ,<{),\|/] (face seal) and [£ x ,e y ,<|>,\|/] 
(cylindrical seal) 

= column vector, transpose of [e] 

= {£} at previous iteration in homing process 

= eccentricity disturbance function (scaler) 

= small increment used in perturbations and in evaluating derivatives 
= angular coordinate, see Fig. 1 
= angular coordinate at start of groove 
= compressibility number, 6pcoR 0 2 /(p 0 C 2 ) 

= groove compressibility number, A3d)a(l-a)sinp 
= viscosity 

= squeeze number, 2AQ 
= global shear stress 
= dimensionless shear stress, xR 0 /(p 0 C) 

= local shear stress 
= rotation about x axis 
= reduced rotation, OR^C 
= rotation about y axis 
= reduced rotation 'FRq/C 
= angular velocity of disturbance 

dimensionless disturbance frequency, Q/co 
= total angular velocity, 0 ), + (o 2 
= dimensionless angular velocity ratio, (o^-coy/co 
= angular velocity of grooved surface 


NAS A/CR— 2004-2 1 3 1 99/VOL3 


xiv 



co 2 = angular velocity of smooth surface 

V = gradient operator, (l/r)(3/30)i + (3/3s)j, on dimensional quantities and 

(l/R)(3/30)i + (3/3 S)j, on dimensionless quantities 
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NOMENCLATURE FOR CODE SPIRALI 


[A] 

[B] 
[B] 

C 

e xA 

e x ,e z 

^xO’^zO 

F 

f 

fa 

fb 

H 

H 

^-BRL 

H, 

^TAP 

h 

h 

h' 


= apparent mass matrix, see Eq. (3-90) 

= damping matrix, see subscripts below for definitions of components 
= dimensionless stiffness matrix, see Eq. (3-89) for component scale factors 
= clearance (cylindrical seal) or reference film thickness (face seal) 

= displacements of center of seal in x,z directions, see Fig. 6 
= dimensionless displacements, e x /C, e/C 
= dimensionless displacement amplitudes, see Eq. (3-60) 

= function defining derivatives, see Eqs. (3-56)-(3-58) 

= shear factor ( l A friction factor) 

= shear factor at surface (a), see Eq. (3-4) 

= shear factor at surface (b), see Eq. (3-5) 

= first-order film thickness 

= dimensionless film thickness, H/C 
= film change defining quadratic "barreling", see Fig. 12 
= first-order film thickness over ridges 
= film change defining linear taper, see Fig. 12 
= film thickness 

= dimensionless film thickness, h/C 
= second-order film thickness 


h' 


h r 

h r 

[I] 

Ic 

If 


= dimensionless film thickness, h'/C 
= film thickness over grooves, see Fig. 10 

= dimensionless film thickness, h/C 
= film thickness over ridges, see Fig. 10 
= dimensionless film thickness, h/C 
= unit diagonal matrix of order N 

= cylindrical seal flag parameter, I c =l for cylindrical seal, I c =0 for face seal 
= face seal flag parameter, Ipl for face seal, lf=0 for cylindrical seal 
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xvn 



Ico 

I 

r,T 

3(Z) 

h 

[K] 

[K] 

[K°] 

[k] 

L 

M x ,M y 

M x ,M y 

m a 

m b 

mo 

N 


n 0 

“e 

P 

P 


P 

P 

p' 

p' 


- rotating groove flag parameter, 1^=1 for grooves on rotor, I 0 =0 for grooves 
on stator 

= unit imaginary number, V-l 
= unit vectors in 0,s directions 
= imaginary part of complex number, Z 
= see Eq. (3-64) 

= stiffness matrix, see subscripts below for definitions of components 

= dimensionless stiffness matrix, see Eq. (3-88) for component definitions 
= stiffness matrix, [K] at E>=0 

= derivative matrix used in numerical solution, see Eq. (3-58) for definitions of 
components 

= seal length, see Fig. 6 
= applied moments about x, y axes 

= dimensionless moments, M^p^), M y /(p 0 To) 

= exponent in shear factor power law relationship for surface(a), f a = n a R™ a 

= exponent in shear factor power law relationship for surface(a), f b = n b R b b 

= exponent for shear factor power law based on Blasius relationship, see Eq. (3-7) 

= order of system of ordinary differential equations, see Eqs. (3-56)-(3-58) 

= number of grooves 

= coefficient in shear factor power law relationship for surface(a), f a = n a R a a 
= coefficient in shear factor power law relationship for surface(a), f b = n b R b b 
= coefficient for shear factor power law based on Blasius relationship, see Eq. (3-7) 
= unit vector normal to grooves, see Fig. 10 
= first-order pressure 
= dimensionless pressure, P/p 0 

= dimensionless first-order local tangential pressure gradient over grooves see 
Eq. (3-94) 

= pressure 

= dimensionless pressure, p/p 0 
= second-order pressure 
= dimensionless pressure, p'/p 0 
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Po 

* 

P 

Pj 

p! 

/V 

Pk 

/V i /V 

Pk>Pk 

Pk»Pk 

Pex 

Pin 

Pex’Pin 

Pg,9’Pg,s 

Pg,0’Pg,s 

Pr,8’Pr,s 

Pr,9’Pr,s 

P 

9 

Q 

^g0’Qgs 

^r8’^lrs 

q n 

q n 

R 

R* 

R„ 

R„ 

r 

r 


= reference pressure 
= dimensionless parameter, pV 0 r 0 /(4C 2 p 0 ) 

= pressure at upstream side of jump in film thickness 
= second-order pressure at upstream side of jump in film thickness 
= component of second-order pressure, p', associated with 5 k , see Eq. (3-65) 

= complex, second-order, dimensionless component pressures, see Eqs. (3-66)-(3-67) 
= complex S dependent factors of p k and p k respectively, see Eq. (3-67) 

= pressure at exit side of seal 
= pressure at inlet side of seal 
= dimensionless pressures, p ex /p 0 , p ir /Po 

= local pressure gradients over grooves in tangential (p g0 ) and transverse (p g s ) 
directions 

= dimensionless local pressure gradients r o p g0 /p o , r 0 p gs /p 0 , see Eqs. (3-24)-(3-25) 

= local pressure gradients over ridges in tangential (p r0 ) and transverse (p r s ) 
directions 

= dimensionless local pressure gradients r o p r0 /p o , r 0 p rs /p 0 , see Eqs. (3-26)-(3-27) 

= power loss 

= dimensionless power loss, ^(Cpo^CD) 

= volumetric flow rate in transverse (s) direction 
= dimensionless flow components over grooves, u g h g , v g h g 
= dimensionless flow components over ridges, u^, v^ 

= flow rate in n p direction relative to spiral groove motion 
= dimensionless flow rate normal to spiral grooves, q n /(h 0 V 0 ) 

= Reynolds number, see Eq. (3-9) 

= characteristic inertia parameter, 2CR/r 0 

= local Reynolds number associated with surface (a), see Eq. (3-6) 

= local Reynolds number associated with surface (b), see Eq. (3-6) 

= radial coordinate, taken as r 0 for cylindrical seal 
= dimensionless coordinate, r/r 0 , taken as 1 for cylindrical seal 
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1*0 


K(Z) 

s 


^in’^ex 

S 


Sj 

s L 


U 

u 


u 

u 


= reference radius, taken as outside radius for face seal and shaft radius for 
cylindrical seal 

= real part of complex number, Z 
= dimensionless coordinate, s/r 0 
= dimensionless coordinate s/r 0 
= dimensionless boundary coordinates s L /r 0 , s R /r 0 
= dimensionless coordinates s in /r 0 , s ex /r 0 

= transverse coordinate, s = r for a face seal and s = z for a cylindrical seal 
= s coordinate at jump in film thickness 

= left boundary coordinate, inside radius for face seal and s L = -L/2 for 
cylindrical seal 

= right boundary coordinate, s R = r 0 for face seal and s R *= L/2 for cylindrical seal 
= s coordinate at exit from seal 
= s coordinate at inlet to seal 
= time 

= dimensionless time, V 0 t/r 0 

= unit vector tangent to grooves, see Fig. 10 
= first order tangential (0) velocity component 
= dimensionless velocity, U/V 0 
= tangential (0) velocity component 
= dimensionless velocity, u/V 0 
= second-order tangential (0) velocity component 
= dimensionless velocity u'/V 0 

= fluid velocity vector 


u. 


u g ,u r 

G g ,G r 


= velocity of surface a, see Fig. 7 

= velocity of surface b, see Fig. 7 
= tangential velocity over grooves (u ), ridges (u r ) 

= dimensionless velocity components Ug/V 0 , u/V 0 
= tangential velocity component at inlet side of seal 
= dimensionless velocity, u in /V 0 
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xx 



A 

U k 

A i A 

u k ,u k 


KA 

v 

V 
V 0 

v 

V 


V 

v' 

V g ,v r 

V g ,v r 


Vj 


Vj 

V{ 

V} 

A 

V k 


v k ,v k 


AA 

w 

T T z 

w 

z 


w x ,w y ,w z 

w x ,w y ,w z 


x,y,z 

{Y} 


| y-new j 


a 

p 


= component of u', associated with 5 k , see Eq. (3-65) 

= complex, second-order, dimensionless component tangential velocities, see Eqs. 
(3-66)-(3-67) 

= complex S dependent factors of uj and u k respectively, see Eq. (3-67) 

= first-order transverse (s) velocity component 

= dimensionless velocity, V/V 0 
= reference velocity 
= transverse (s) velocity component 
= dimensionless velocity component, v/V 0 
= second-order transverse (s) velocity component 
= dimensionless velocity v'/V 0 
= transverse velocity over grooves (v ), ridges (v r ) 

= dimensionless velocity components v/V 0 , v/V 0 
= s velocity component at upstream side of jump in film thickness 
= dimensionless velocity v/V 0 

= second-order s velocity component at upstream side of jump in film thickness 
= dimensionless velocity v]/V 0 
= component of v', associated with 6 k , see Eq. (3-65) 

= complex, second-order, dimensionless component transverse velocities, see Eqs. 
(3-66M3-67) 

= complex S dependent factors of v k and v k respectively, see Eq. (3-67) 

= first-order axial load applied to face seal 

= dimensionless load, 

= applied loads in x, y, z directions 
= dimensionless loads, w^p^), w^p^), w z /(p 0 To) 

= coordinate variables, see Fig. 6 

= column vector of independent variables at old grid point, see Eqs. (3-56)-(3-58) 
= column vector of independent variables at new grid point, see Eqs. (3-56)-(3-58) 

= groove to pitch ratio, A0g/(A0 g +A0 r ) 

= spiral groove angle, angle between grooves and surface velocity 
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Ah 

Ah 

Ap 

Ap 

Ap' 

Ap g ,Ap r 

Ap g ,Ap r 

AS 

A0 

A0 r 

8 

5 

S k 


e k 

c 

T1 k 

0 

X 

p 

p 

X 

x 
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- jump in film thickness (downstream - upstream) 

= dimensionless jump in film thickness, Ah/C 

= jump in pressure at s j5 p-p ; 

= dimensionless jump in pressure, Ap/p 0 

= dimensionless second-order jump in pressure at Sj, (p'-pJVPo 

= spiral groove downstream - upstream pressure jump at entrance to grooves, 
(Ap g ) or ridges, (Ap r ) 

= dimensionless spiral groove pressure jumps, Apg/p 0 , Ap/p 0 

= variable increment between grid points used in numerical solutions, see 
Eqs. (3-57)-(3-58) 

= circumferential extent of groove-ridge pair, see Fig. 10 
= circumferential extent of groove, see Fig. 10 
= circumferential extent of ridge, see Fig. 10 
= groove depth, h g -h r 

= dimensionless groove depth, 8/C 

= dimensionless component film displacements, see Eqs. (3-59)-(3-61) 

= complex dimensionless component displacements associated with and -Q, 
respectively, see Eq. (3-63) 

= set of amplitudes, see Eq. (3-64) 

= loss coefficient due to contraction 

= see Eq. (3-64) 

= angular coordinate, see Fig. 6 

= correction factor applied to x c , see Eq. (3-94) 

= viscosity 

= loss coefficient, see Eq. (3-18) 

= density 

= tangential shear stress 
= dimensionless shear stress, r 0 x/(Cp 0 ) 

= tangential shear stress, surface (a) 

= dimensionless shear stress, r 0 x/(Cp 0 ) 
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= shear traction vector for surface (a), see Fig. 7 

= tangential shear stress, surface (b) 

= dimensionless shear stress, ^^(Cpo) 

= shear traction vector for surface (b), see Fig. 7 

= dimensionless tangential Couette shear stress, (x a + x b )/2 

= Couette shear traction vector, (x a + x b )/2 

= dimensionless tangential Poiseuille shear stress, (x a - x b )/2 

= Poiseuille shear traction vector, (x a - x b )/2 

= dimensionless effective shear stress, see Eq. (3-92) 

= values of x* obtained from first-order solutions for grooves and ridges, 
respectively 

= turbulent shear functions, see Eqs. (3-14)-(3-15) 

= global functions defined by Eqs. (3-42)-(3-43) 

= turbulent shear functions over grooves, see Eq. (3-22) 

= turbulent shear functions over ridges, see Eq. (3-23) 

= function for predicting pressure change at jump in film thickness, see Eq. (3- 
= rotation about y axis, see Fig. 6 

= dimensionless amplitude, see Eq. (3-60) 

= angular velocity of oscillation 

= dimensionless angular velocity, Qr/V 0 

= see Eq. (3-71) 

= angular velocity of rotor 
= dimensionless angular velocity, cor/Vo 


Subscripts 

a,b = relating to surface (a),(b) shown in Fig. 2, generally (a) is moving, (b) is 

stationary 

ex = exit 

g = grooves 
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H,U,V 


partial derivatives with respect to H, U and V, respectively 
inlet 


in 

J = jump discontinuity in global film thickness 

k = index associated with displacements, see Eq. (3-64) 

L = left or lower value of s 

R = right or upper value of s 

r = ridges 

x = pertains to force or displacement in the x direction 

y = pertains to force or displacement in the y direction 

§ = pertains to moment or rotation about x axis 

\|/ = pertains to moment or rotation about y axis 
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1-0 INTRODUCTION 


NASA's advanced engine programs are aimed at progressively higher efficiencies, greater 
reliability, and longer life. Recent studies have indicated that significant engine performance 
advantages can be achieved by employing advanced seals [1]*, and dramatic life extensions 
can also be achieved. Advanced seals are not only required to control leakage, but are 
necessary to control lubricant and coolant flow, prevent entrance of contamination, inhibit the 
mixture of incompatible fluids, and assist in the control of rotor response. 

Recognizing the importance and need of advanced seals, NASA, in 1990, embarked on a five- 
year program (Contract NAS3-25644) to provide the U.S. aerospace industry with computer codes 
that would facilitate configuration selection and the design and application of advanced seals. 

The program included four principal activities: 

1 . Development of a scientific code called SCISEAL, which is a Computational Fluid 
Dynamics (CFD) code capable of producing full three-dimensional flow field 
information for a variety of cylindrical configurations. The code is used to enhance 
understanding of flow phenomena and mechanisms, to predict performance of complex 
situations, and to furnish accuracy standards for the industrial codes. The SCISEAL 
code also has the unique capability to produce stiffness and damping coefficients that 
are necessary for rotordynamic computations. 

2. Generation of industrial codes for expeditious analysis, design, and optimization of 
turbomachinery seals. The industrial codes consist of a series of separate stand-alone 
codes that were integrated by a Knowledge-Based System (KBS). 

3. Production of a KBS that couples the industrial codes with a user friendly Graphical 
User Interface (GUI) that can in the future be integrated with an expert system to 
assist in seal selection and data interpretation and provide design guidance. 

4. Technology transfer via four multiday workshops at NASA facilities where the results 
of the program were presented and information exchanged among suppliers and users 
of advanced seals. A Peer Panel also met at the workshops to provide guidance and 
suggestions to the program. 

This final report has been divided into separate volumes, as follows: 

Volume 1: Executive Summary and Description of Knowledge-Based System 

Volume 2: Description of Gas Seal Codes GCYLT and GFACE 

Volume 3: Description of Spiral-Groove Codes SPIRALG and SPIRALI 

Volume 4: Description of Incompressible Seal Codes ICYL and IF ACE 

Volume 5: Description of Seal Dynamics Code DYSEAL and Labyrinth Seal Code KTK 

Volume 6: Description of Scientific CFD Code SCISEAL. 


*Numbers in brackets indicate references located in Section 4.0. 
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This volume describes two industrial codes used to predict the behavior of both cylindrical and 
face seals with or without the inclusion of spiral grooves: SPIRALG and SPIRALI. This includes 
load capacity (for face seals), leakage flow, power requirements and dynamic characteristics in the 
form of stiffness, damping, and apparent mass coefficients in 4 degrees of freedom for cylindrical 
seals and 3 degrees of freedom for face seals. The code SPIRALG treats laminar flow in gas 
seals operating at finite eccentricities. The code SPIRALI deals with high-speed incompressible 
flow and includes turbulence and other inertia effects but is limited to small eccentricities. 
References 2 and 3 provide the details of code implementation. 

Turbulence and inertia have been included in the incompressible code, SPRIALI, at the 
expense of the effects of large eccentricity. Cryogenic liquids have relatively low kinematic 
viscosities, and seals for such liquids are sometimes designed to run at relatively high 
clearances. In such cases, stiffness effects associated with inertia can provide far greater 
stabilization than eccentricity effects for laminar gas seals. In addition, SPIRALG can be run 
at high ambient pressure to simulate noncavitating incompressible seals at finite eccentricities. 
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2.0 INDUSTRIAL CODE SPIRALG - GAS-LUBRICATED, 
SPIRAL GROOVE, CYLINDRICAL AND FACE SEALS 

The computer code SPIRALG predicts performance characteristics of gas-lubricated, spiral 
groove, cylindrical and face seals. Performance characteristics include load capacity, leakage 
flow, power requirements and dynamic characteristics in the form of stiffness and damping 
coefficients in four degrees of freedom for cylindrical seals and three degrees of freedom for 
face seals. These performance characteristics are computed as functions of seal and groove 
geometry, loads or film thicknesses, running speed, fluid viscosity, and boundary pressures. 

The basic assumptions that have gone into the computer code are listed below: 

1. The flow is assumed to be laminar and isothermal. 

2. Inertial effects are neglected. 

3. The gas is assumed to be ideal. 

4. The film thickness is assumed to be small compared with seal lengths and diameters 
but large compared with surface roughness and the mean free path of the gas. 

5. Narrow groove theory is used which characterizes the effects of grooves by a global 
pressure distribution without requiring computations on a groove by groove basis. 

This involves neglecting edge effects and local compressibility effects associated with 
groove to groove pressure variations. In general, narrow groove theory is valid when 
there are a sufficiently large number of grooves so that 2n/N g « 1, where p is the 
groove angle and N g is the number of grooves. 

6. Transient effects are treated with the use of small perturbations on a primary steady- 
state flow. These transient effects are characterized by stiffness and damping 
coefficients that are dependent on the disturbance frequencies. 

7. Although displacements and misalignments are treated, machined surfaces for face 
seals are assumed to be flat and machined clearances for cylindrical seals are assumed 
to be constant. 

The above assumptions still leave the code applicable to a broad range of applications. Seals 
generally have small clearances and gasses have low densities resulting in sufficiently low 
Reynolds numbers for laminar flow. Practical designs should contain a fairly large number of 
grooves to ensure smooth, isotropic operation. At high sealed pressure differences the flow 
could become sonic thereby invalidating the first two assumptions but this will usually n'dt f>e 
the case and can readily be checked based on the predicted leakage flow. Elastic and thermal 
distortions as well as machining tolerances should also be estimated to validate the constant 
clearance assumption. The overall accuracy of the program will depend on the grid size used. 
Factors such as high compressibility or squeeze numbers, small values of the minimum film 
thickness to clearance ratio and large values of the length to diameter ratio could require 
either a large number of grid points or carefully selected variable grids. 
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2.1 Theoretical Development 

The first formulation of the equations governing gas-lubricated spiral groove bearings is generally 
credited to Vohr and Pan [4]. A more concise formulation is given in a second report [5] by 
these authors that has been used by Smalley [6] as a starting point in his generalized numerical 
treatment of the performance of spiral groove gas bearings. The work performed by Smalley may 
be applied to both bearings and seals. A principal limitation in all of the above references relates 
to the fact that solutions have only been provided for one dimensional forms of the equations 
which have been obtained by linearizing them based on near concentric and aligned conditions. 
The work described here deals with the numerical solution of the nonlinear equations for gas- 
lubricated spiral groove seals at both eccentric and misaligned conditions. 


2.1.1 Formulation of Equations Governing Gas-Lubricated Spiral Groove Seals 

For completeness, a derivation of the narrow groove equations for spiral groove gas bearings 
and seals along the lines of that developed in Reference 5 will be provided here. Coordinate 
variables will be used to make the equations applicable to both cylindrical and face seals as 
can be seen with the aid of Figure 1. The circumferential coordinate, 0, is as shown in 
Figure 1. The transverse coordinate is described by the variable, s, which is taken to equal 
the radial coordinate, r, for a face seal and the axial coordinate, z, for a cylindrical seal. The 
quantity r, when it appears will denote radial position for a face seal and should be set equal 
to the shaft radius, Rq for a cylindrical seal. 

The isothermal, compressible form of the "Reynolds" equation may be written as a flow 
balance equating the divergence of the flow vector, q', to the flow per unit area squeezed out 
by the time rate of decrease of the film thickness, 


^■q' = — — (rq.. 7 ) + = 

r as r 00 


Qa 


2-1 


The local flow vector q' = qji + q'j represents the mass flow rate per unit transverse length 
divided by the density at a reference pressure, p 0 , which may be written in vector form as 


q 7 


Po 2 Po 
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Since surface motion will be in the circumferential direction, the surface velocity vectors may 
be written as u, = rco,i and u 2 = r(i) 2 i and the components of q' become 


h 3 1_ + r (Qj + 0) 2 p 7 h 

12|i p Q r 00 2 p 0 


2-2 


h 3 p ; 0p ; 

12m Po & 
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The "squeeze film" term or displaced mass flow per unit area due to film motion, divided by 
the density at p 0 is 


q 


/ 

A ~ 


_ 1 flp'h) 

Po 31 


2-4 


One could substitute Equations (2-2) - (2-4) for the corresponding flow quantities in Equation 
(2-1) to obtain the usual form of the compressible Reynolds Equation which could in principle 
be solved, for any film thickness profile, h(s,0) and appropriate boundary conditions, for the 
pressures or flow components to obtain the pressure distribution. These could in turn be 
integrated to obtain the various forces and moments associated with the given bearing 
geometry. The torque opposing the motion of say the smooth surface may be determined, 
once the pressure distribution is known, by integrating the shear stress relationship that arises 
in the development of Reynolds equation 


/ = h 0p 7 

2r 00 



-<*i 

h 
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The difficulty encountered in obtaining full numerical solutions to the above equations relates to 
the complexity of the grid network necessary to adequately describe the geometry of a surface 
containing the large number of spiral grooves usually required to provide sufficiently smooth 
pressure distributions to make the load characteristics independent of whether shaft displacement 
is over a ridge or over a groove. Narrow groove theory is generally used to circumvent this 
difficulty (References 4-6). It will be implemented here, as well and is described below. 

Narrow groove theory provides the limiting form of the solution to Equations (2-1) - (2-5) as 
the number of grooves, N g , becomes large, with the groove angle, p and the groove to pitch 
ratio, a, held constant. The discontinuities in film thickness associated with the grooves will 
give rise to discontinuities in the pressure gradients at the ridge-groove interfaces as illustrated 
schematically in Figure 2. The local pressure profile p' is shown by the sawtooth lines 
whose lower vertices, for purposes of illustration, are connected by the "global" pressure 
profile, p. The global pressure profile does not necessarily lie at the lower vertices of the 
local pressure profile but could lie anywhere between the lower and upper vertices. In the 
limit as the number of grooves becomes large the curve connecting the upper vertices will 
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Figure 2. Schematic of Spiral Groove Parameters, Global and Local Pressures 
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approach the curve connecting the lower vertices. This limiting behavior is not true of 
3p730, 3p73s or h, which will have different values over lands and grooves no matter how 
large the number of grooves. Narrow groove theory requires the development of expressions 
for the local (primed) quantities in terms of global quantities that approach single valued 
limits as the number of grooves becomes large. The local film thickness and pressure 
derivatives over the grooves will be denoted by h g , (3p730) g and (3p73s) g respectively and 
by 1\, (3p730) r and (3p73s) r over the ridges. (The subscript r has been used here to denote 
ridges for consistency with References 4-6 and should not be confused when used in a 

different context later in the report to denote the right-hand boundary pressure or with the 

radial position variable, r, which is not used as a subscript.) 

When the number of grooves becomes large, the sawtooth portion of the local pressure 
variation may be approximated with linear representations as shown in Figure 2. Thus, 
equating pressures over a groove-ridge pair in the circumferential direction 

a p _ Ap g ' ^ Ap r ' (ap'i A0 g fapO A 0 r 

A0 A0 A0 ia0j g A0 iaeJ r A0 

Noting that A0 g /A0 = a and A0/A0 = 1 - a and replacing Ap/A0 with 3p/30 as A0 — » 0, the 
above equation becomes 



+ (1 -a) 



The corresponding relationship in the transverse direction, 


2-6 



is obtained in a similar manner. 

The remaining two equations required to solve for the four local pressure derivatives are 
obtained from continuity considerations. First, the pressure must be continuous at each 
groove-ridge interface, thus the derivative of the pressure in the direction of the interface, 
Vp'«tp, must also be continuous. The second requirement is for continuity of the flow 
normal to each groove-ridge interface as measured in a frame of reference moving with the 
grooves, (q' - rco 1 hp7p 0 i)*n p . The unit tangent and normal vectors for a logarithmic spiral 
are given by 

t p = cospf + sinpj , n p = sin(3f - cospj . 

The first of the above conditions requires continuity of 


+ sinp^ 
r 50 0s 
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The second condition requires continuity of 
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One may substitute Equations (2-2) and (2-3) for the circumferential and transverse 
components of the flow vector at each groove-ridge interface, respectively to obtain 


3 r 
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Density variation, p7p 0 , is continuous at each interface and cancels out of Equation (2-9). 

Equations (2-6) - (2-9) represent the four linear equations needed to solve for the local 
pressure derivatives. We may obtain the solution by first solving Equations (2-8) and (2-9) 
for the components of the local pressure gradient over the grooves in terms of those over the 
ridges. The resulting equations may be written as 


1 _ 

r 
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2-11 
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One may now substitute Equation (2-10) for (8p730) g in Equation (2-6) and Equation (2-11) 
for (3p73s) g in Equation (2-7) to obtain two linear equations for the components of the ridge 
pressure gradient which may in turn be solved to yield the following expressions: 


d?!) _ 

dQ 


[h a 3 -«(h a 3 -h r 3 cos 2 p)]i^ - a(h 3 -h 3 )sinpcosp|! - 6^r( 

r do ds 

(1 - a)h 3 + ah r 3 


« 2 - co 1 )a(h - h r )sin 2 p 
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(1 - a)h 3 * ah 3 
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The components of the local groove pressure gradient may be expressed in terms of the above 
ridge components by simple rearrangement of Equations (2-6) and (2-7): 


l 3p, l 
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Now that expressions have been developed for the components of the local pressure gradients 
in terms of global ones, it is necessary to determine the global flow components q s and q 0 
and a global squeeze film term q A that may be substituted for the local ones in the flow 
balance given by Equation (2-1). These global flow components are determined by matching 
mass flow rates over a groove-ridge pair with the mass flows obtained by integration of the 
local flow components over the same interval. 

If 0 g is taken as the circumferential coordinate at the start of a groove, the transverse flow 
crossing an arc at fixed s, subtending a groove-ridge pair in the interval 0 g < 0 < 0 g +A0 is 
given by the left-hand term in the relationship 
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The approximation to the integral in the above expression was obtained by dividing the integration 
interval, A0 into sub-intervals for the groove, A0 g and ridge, A0 r and approximating q' s , noting 
that as the number of grooves becomes large 0p73s, will approach a constant value within each 
sub-interval. Since the pressure at the groove-ridge interface is continuous, the local density 
variation term, p'/p 0 was replaced by its global value p/p 0 . The far right-hand term in the above 
expression is based on the definition of the transverse component of the global flow rate described 
above. The right two equalities may be solved for q s as 


q s 


12p Po 
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One may obtain a relationship for the circumferential flow component q e in a similar manner 
by integrating q' e , given by Equation (2-3), at fixed 0, over a groove-ridge pair 
( s g < s < s g +As ), approximating the integral over each sub-interval as above and equating the 
result to q e As. The resulting expression may be written as 


q 9 
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By integrating the squeeze film term q/, given by Equation (2-4), over an area rA0As, 
equating it to (^^©As and noting that the groove area fraction will be a and the ridge area 
fraction will be 1 - a, the following expression is obtained: 

q* = !(p[«V(i -«)hj) . 2-18 

Po at ' 


The global shear stress x, may be determined by integrating the local shear stress x', given by 
Equation (2-5), with respect to 0 over the interval 0 g < 0 < 0 g +A0, invoking the narrow 
groove approximations and equating the result to xA0. The resulting expression may be 
written in the form 
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Equation (2-1) may now be applied directly to the global flow vector q = q e i + qj, as 
V«q = and by substituting Equation (2-18) for q A and putting the result in dimensionless 
form one obtains: 


VQ 
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R 00 
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The components of the global flow vector, q 0 and q s are given in terms of the local pressure 
derivatives by Equations (2-16) and (2-17), respectively. These local derivatives are, in turn, 
given in terms of the global ones by Equations (2-12) - (2-15). The global flow components 
may be expressed completely in terms of global pressure derivatives by first substituting 
Equations (2-14) and (2-15) for the local pressure derivatives over the grooves and then 
substituting Equations (2-12) and (2-13) for the local pressure derivatives over the ridges. 
One may then collect terms and put the resulting two equations in dimensionless form to 
obtain the following expressions for the components of the dimensionless flow vector 

Q = Qe* + Qj: 
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Q s - -(1 +P) 
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The dimensionless variables associated with the above equations are 
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The dimensionless gage pressure P in the above equations is taken relative to the absolute 
pressure p 0 which will henceforth be taken as the minimum of the two boundary pressures in 
absolute units. The dimensionless parameters associated with the above equations are 
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and the column matrix containing spiral groove coefficients, kj(ot,p,r), in the above equations is 


k = 


«(i -«)(r 3 -ifsin 2 p +r 3 
(1 - a)r 3 + a 

a(1 -a)(r 3 -1) 2 sinpcosp 
(1 - a)r 3 + a 

g(1 -a)(r 3 -1) 2 cos 2 p +r 3 
(1 - a)r 3 + a 

( r 3 -i) 

(1 - a)r 3 + a 

(1 - a)r + a 
f 


(T -1)sinp 
(1 - a)r 3 + a 

«(1 - a)(r 3 - 1)(r -1)sinpcosp 
(i - «)r 3 + « 

«(i -q)(r 3 -i)(r-i)cos 2 p +«r + (i -g)r 3 
~ (1 - a)r 3 + a 
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Only the first four components of k are used in Equations (2-21) - (2-22). The remaining 
components are used in evaluating the shear stress. The relationships for k„ i=l,2,3>4 derived 
here are consistent with Equation (3.27) of Reference 5. 

The global shear stress is obtained by substituting Equation (2-14) for (0p'/30) g /r in Equation 
(2-19) then substituting Equation (2-12) for (3p730)/r in the resulting expression. The latter 
result may be expressed in dimensionless form as 


A toR 

^TT 
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k .^1 

^ r ae 


2-26 


which is consistent with Equation (3.88) of Reference 5. 

The equations presented thus far are directly applicable to either a cylindrical seal or a face seal. 
As mentioned earlier, a face seal is represented in the above equations by setting the transverse 
coordinate s equal to the radial coordinate r. This is equivalent to setting S = R in dimensionless 
form. A cylindrical seal is represented in dimensionless form by setting S = Z and R = 1. 
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The quantities required to characterize the groove dimensions are shown in Figure 2. If by 
convention co is taken to be positive (surface motion in the direction of increasing 0), then the 
groove angle, p, will be the angle measured from the groove to the direction of surface 
motion associated with to. A positive acute value of p will tend to pump in the positive S 
direction if the grooves are on the stator and in the negative S direction for grooves on the 
rotor. By setting the groove depth parameter 8 = 0, T becomes 1 and Equations (2-20) - (2- 
26) reduce to those for ungrooved seals. By treating k and 8 as sectionally continuous 
functions of S, these equations may be applied to composite smooth and grooved geometries 
with Q s and P held continuous at all transition boundaries. 

The film thickness relationship for H r , which may be applied to either a cylindrical seal or a 
face seal is 


H r = 1 - € z - (e x +i|tS)cos0 - (e y - <J>S ) sin 0 2-27 

with s z = 0 for a cylindrical seal and e x = e y = 0 for a face seal. 

The boundary pressures will be taken to be p, and p r at the inside and outside radii 
respectively for a face seal or at the two ends (z = -L/2 and Z = L/2) for a cylindrical seal. 
This is expressed in dimensionless form 

P = P, at S = S, , P = P r at S = S r . 2-28 

The remaining boundary condition relates to periodicity with respect to 0 which requires P 
and Q 0 to have the same values at 0 = 0 as they do at 0 = 2k: 

P la=o = P la-2, and Q 8 U =Q ele-2, • 2 ' 29 

The above treatment is intended to represent a complete statement of the mathematical 
problem for determining the pressures and surface shear stresses in plain or spiral groove face 
or cylindrical seals. The rest of this section will deal with the numerical determination of the 
pressure distribution and the computation of related quantities such as loads, leakage, power 
loss, stiffness and damping. 


2.1.2 Discretization of Pressure Equations 

Discretization will be carried out with the use of the cell method [7] which involves the 
performance of a flow balance about each interior grid point. One may integrate Equation (2-20) 
over an arbitrary control area within a seal 

f V-QdA + f 4 [(«S + H r )(1 + P)]dA = 0 
J A J A at 
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and apply the divergence theorem to the first integral on the left to obtain the relationship 


<£ Q-ndS + f — [(<x6 +H r )(1 +P)]dA = 0 , 2-30 

J S J A a 

which will be used as a starting point in the discretization process. 

A grid network may be set up along with flow control areas about each grid point as shown in 
Figure 3. The grid will contain M lines in the S direction including boundaries and N lines in the 
0 direction from 0 = 0 to 0 = 2 tc, inclusive. The grid points at the intersections of these lines are 
noted by the solid circles. Flow control areas to be used in evaluation of the integrals in 
Equation (2-30) are set up about each grid point as shown by the shaded area in Figure 3. The 
comers of the flow control area denoted by the shaded points marked 1,2,3, and 4 are located at 
the geometric centers of the rectangles formed by the grid lines and will be referred to as half 
grid points. The flow components labeled Q| 2 etc., represent the components of the flow vector 
in the positive coordinate directions as indicated by the arrows. The subscripts (12 etc.) refer to 
the line connecting points 1 and 2, and the superscripts (+,-) refer to the positive or negative side 
of the point of intersection with the grid line. 

We will adopt the convention that the subscripts i,j refer to grid points and subscripts such as 
i+ 1 / 2 ,j + 1 /2 refer to half grid points. The value of the radius R at half grid point 2 would thus 
be R i+I/j . The differential control length, dS, in Equation (2-30) will be approximated by the 
lengths of the various lines or arcs bounding the flow control area thus AS| 2 refers to the 
length of the line associated with Q} 2 described above which for this example would be AS/2. 
Similarly the arc length associated with Qj 4 would be AS U = Ri^AO^/2. 

The area element dA will be made up of the parts of the flow control area in each of the four 
quadrants about the center point (i,j) and numbered based on the shaded half grid point that 
each contains. Thus AAj = (Rj+R^AOjAS/S, etc. and the discretized form of Equation (2-30) 
may be written as:The flow components on the left-hand side of Equation (2-31) are obtained 

Qi 2 AS 12 + Q 12 AS 12 + Q 14 AS 14 + Q 14 AS 14 -Q 34 AS 34 -Q 34 AS 34 -Q 23 AS 23 -Q 23 AS 23 = 
at v 

2-3 


(a6 +H r ) i4 , jti AA, +(<x5 +H r ) uij _ 1 AA 4 +(a6 + H r ) i _ 1 AA 3 +(«6 +H r ) jJ . j)d AA 2 

2 2 2 2 2 ,J 2 2 1 2 J 


from discretization of Equations (2-21) and (2-22). A numbering system for the 9 pressures 
at the grid point (i,j) and the 8 surrounding points is shown in Figure 3, where P 1 = P i+1 j_j, 

P 5 = Pjj etc. The determination of the flows out of the sub-area containing the half grid point 
labeled 1 is discussed here as an example. The flow component Q} 2 , is determined from 
Equation (2-21). The derivative of the pressure normal to the line connecting points 1 and 2 
is evaluated at the 
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intersection of the line with the grid line. The tangential derivative is evaluated at the half 
grid point, 1 as the average of the difference between P 3 and P 6 with that between P 2 and P 5 , 
divided by AS;. Thus, 


1 9P _ p 6 - p 5 

R 30 R| A0, 


9P 

as 


( p ,-p,)*( p 2 - p 5 ) 

2AS| 


(for q; 2 ) . 


The flow component Qj 4 is determined in a similar manner from Equation (2-22). The 
normal derivative is approximated as the difference between P 2 and P 5 divided by AS; and the 
tangential derivative is approximated as the average of the differences between P 3 and P 2 , and 
P 6 and P 5 , divided by R i+1/2 A0j. 


3P B P 2 — P 5 

as ASj 


j_ap 

r ae 


( P 3 — P 2 ) "*" ( P 6 — P 5 ) 

2R,.,A0 | 

2 1 


(for q; 4 ) . 


For both of the above flow components the pressure in the (1+P) term appearing in Equations 
(2-21) and (2-22) is evaluated at the half grid point by averaging the four surrounding 
pressures (P 2 +P 3 +P 5 +P 6 )/4. All of the remaining quantities (R, H r , k lf k 2 , k 3 , k 4 , a, (3 and 8) 
are evaluated directly at the half grid point. 

The flow balances over the other quadrants are performed in a similar manner. For steady- 
state conditions, the right-hand side of Equation (2-31) will be 0 and the flow balance about 
any interior grid point (i,j) may be written in the form 


F ij ( H r i P 1 i P 2 > P 3 > P 4 ’ P 5 ’ P 6 ’ P 7 ’ P 8 ’ P 9 ) = 0 


2-32 


The definition of Fy may be extended to make Equation (2-32) applicable to the ends of the 
seal as well as the interior points by applying Equation (2-28) at the endpoints as follows: 


p ii = p 5 - P i (i-1) and F Mi =P 5 -P r (i = M) 


The solution to Equation (2-32) may be used to provide all of the steady-state quantities such 
as pressures, forces, moments, flow rate and power loss. The inclusion of the right-hand side 
of Equation (2-31) will be necessary for determination of frequency dependent stiffness and 
damping coefficients which will be discussed later. 


2.1.3 Newton-Raphson Linearization Procedure 

The Newton-Raphson [8] procedure is perhaps the most widely used method for obtaining 
solutions to nonlinear systems of algebraic equations and is described in many textbooks on 
numerical methods such as Reference 8. A procedure similar to that used here is described in 
a paper by Artiles, Walowit and Shapiro [9]. 
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The procedure is started with an initial pressure distribution that satisfies the end conditions 
given by Equation (2-28). A new set of approximations to the pressures in Equation (2-32), 
P k ncw may be obtained by linearizing Fy about a previously established set of approximations 
P k as follows: 


F„ 



= 0 


2-33 


where a forward difference 


dp k 


Fy ( H ( , P, , P k +t) 


■~.P 9 ) ~ F, j (H r ,P 1 P 9 ) 


may be used to numerically evaluate the partial derivatives. 

Pressures without the superscript "new" relate to the previous or "old" approximation. It 
should be noted that the function Fy will not be 0 unless the pressures comprising its 
arguments are exact. If we go back to using grid notation for P (P ;j in place of P 5 etc.) and 
introduce the column vector {Pj new } as the M new pressures at the jth column of grid points, 
Equation (2-33) may be written in the following form: 


[C l ]{P i new } + [E'KPjT) + [ D J ] { P,"i ew } = { R ' } , 


2-34 


where [C j ], [E j ] and [I>] are tri-diagonal matrices whose interior elements, from Equation (2- 
33), are 


C j 


0F ;; 


i,i+k 


0P 


, E 


0F :i 


i,i+k 


i+k,j 


0P 


, D 


3F :i 


i,i+k 


i +k,j -1 


3P: 


, k = -1,0,1 ; i = 2,..., M -1 . 


i +k,j +1 


The interior elements of the column vector {R J } are 


R i' = E ( C M* P i*,i + E M* P i*,H + - F i, 


k =-1 


The above equations may also be applied to the comer elements to produce the result 


C 


1,1 


- r * 


- 1 


Eli - E* iM - - Dm ; m _ 0 i — Pj , Ri - P r 


Equation (2-34) represents a linear system of simultaneous equations that may be solved by 
various matrix inversion procedures. The method used here is the column or transfer matrix 
method, which is described in References 7 and 10. It has been used extensively in solving 
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finite difference problems associated with various forms of the lubrication equations and 
produces accurate results in a fairly efficient manner. Convergence of the Newton-Raphson 
procedure is generally obtained within three to six iterations, depending on degree of 
nonlinearity and the accuracy required. 


2.1.4 Determination of Loads, Moments, Torque, and Leakage 

The dimensionless loads and moments may be obtained by integrating the pressure 
distribution over the seal area as shown below: 


2 s S r 2s S r 2s S r 

W x = f f PcosSRdSde, W y = f f Psin0RdSd0, W 2 = J J PRdSd0, 

OS, OS, OS, 

2-3 

2s S r 2s S r 

M x = -J J Psin0RSdSd0 , M y = J J Pcos0RSdSd0 . 

OS, OS, 

The dimensionless torque is obtained from integration of the shear stress given by Equation 
(2-26) over the seal area: 


2s s r 

f =sign(Q) J | xRdSd0 . 

0 s, 


2-36 


The sign(cb) term has been added to make the torque positive when it opposes the net surface 
motion regardless of which surface (smooth or grooved) is moving. 

Finally, the dimensionless leakage flow, Q in going into the seal at S = S! may be obtained 
from integration of Q 0 , given by Equation (2-21), over the circumference of the seal: 


2s 

Q, n = J Q„Rd0 . 

0 


2-37 


The integrand in the above expression is evaluated by summing the flow components to the 
right of the first 0 grid line in the same manner as that used in developing Equation (2-31). 
It should be noted that any value of S can be used since Q in is independent of S. 

The physical quantities corresponding to the dimensionless ones given above are 
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W{x,y,z} Rq Po^{x,y,z} 


M 


{x,y} 


R 0 p 0 M 


{x.y} 


q ‘ n= i^ Qin ’ T=R » p ° ct 


2-38 


In the above equations, the loads W x and W y apply only to a cylindrical seal and W z applies 
only to a face seal. The leakage flow q in is the volumetric flow rate going into the seal 
measured at pressure p 0 . 


2.1.5 Determination of Stiffness and Damping Coefficients 

Equation (2-20) with flow components Q e and Q s given by Equations (2-21) and (2-22) 
represents a second-order nonlinear partial differential equation that may be used to define a 
second-order nonlinear operator G, such that 

G(P,H r ) = -—[(<*5 +H r )(1 + P)] . 2-39 

at 


The determination of P under steady-state conditions, where the right-hand side of Equation 
(2-39) is 0, was described earlier in this section. These steady-state pressures will now be 
referred to as P. The various eccentricities and rotations used in determining H r from 
Equation (2-27) may, for convenience, be put in the form of a row matrix as 


[e z ,(|>,<|r] , 
[e x ,e y ,4>,<|j] 


(face seal) 
(shaft seal) 


2-40 


and Equation (2-27) may be written as 

H r = 1 +[e]{a} 

where the column vector {a} is given by 


{a} 


{ -1 ,Ssin0, -Scos0} , (face seal) 

{ -cosQ, -sin0,Ssin0, -Scos0} , (shaft seal) 


2-42 


One could develop a perturbation analysis for prediction of stiffness and damping coefficients 
with the following procedure: (a) perturb say the ith component of the eccentricity matrix in 
Equation (2-40) by qe', where q is a small parameter and e' is time dependent; (b) express 
H r in the form H r = H + qe'fa} and the corresponding pressures as P = P + q{P' }; (c) 
substitute the above expressions for P and H r in Equation (2-39); (d) expand the resulting 
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expression neglecting terms of order q 2 and higher; (e) collect terms of order r\. The resulting 
expression could be written in the form: 


S£{P'( + {b}e' = -(«5 + H)2ilp. - (1 + P){a}-^ , 

at at 


where ££ is a second-order linear operator given by 


C£ = a, + — + Ag-— — + A 4 -A + A 6 A + A 6 

'as 2 asae 3 ae 2 as 5 ae 6 


2-43 


2-44 


The coefficients in the above equations (A,, A 2 , {b} etc.) will depend on the coordinate variables 
as well as P, H and their various derivatives. Only the form of the above equations is important 
to the numerical procedure under development and the significant amount of algebraic 
manipulation required to determine these coefficients will be shown to be unnecessary. 

If the time dependence of the eccentricity is restricted to oscillatory disturbances one may set 
e' = e 3m and look for solutions in the form { P' } = {P*}e 3ot , where {P*} is complex but in- 
dependent of time. When this transformation is introduced into Equation (2-43), the result is 

SEJP’l + { b} = -go[(aS + H){P *} + (1 + P){a}] . 2 ' 45 

The representation of {P*} and the eccentricity coefficients {a} as column vectors relates to 
the fact that each of the eccentricities must be perturbed to obtain the compete stiffness matrix 
but Equation (2-45) is solved independently for each perturbation. The periodic boundary 
conditions given by Equation (2-29) also apply to Equation (2-39) (continuity of {P*} and 
3{P*}/d0 is sufficient when H r and the spiral groove coefficients k 1? ...,k 4 are continuous 
functions of 0 as they are here). The end boundary conditions given by Equation (2-28) 
become {P*} = 0 at S = S, and S = S r . 

If Equation (2-45) were solved for {P*} subject to the above boundary conditions, all of the 
dimensionless stiffness and damping coefficients could be obtained by substituting {P*} for P 
in Equation (2-35). The real parts of the computed forces and moments would be in phase 
with the eccentricity perturbations and constitute the dimensionless stiffness coefficients. 

Thus K yx would correspond to the real part of W y computed from the component of {P*} 
associated with the perturbation in e x and K^ y would correspond to the real part of M x 
computed from the component of {P*} associated with the perturbation in e y etc. In a similar 
manner, the dimensionless damping coefficients which are 90° out of phase with the 
eccentricity perturbations would be obtained by dividing the imaginary parts of the forces and 
moments computed in the manner described above by a. 

The parameter a, is a dimensionless disturbance frequency referred to as the "squeeze 
number" and is given by o = 2AQ/0) where Q is the angular velocity of the disturbance. The 
limiting form of the stiffness and damping coefficients as a — » 0, is of interest as it applies to 
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incompressible flow, and the limiting stiffnesses are used in the homing procedure that has 
been implemented for determining eccentricities from given loads which will be described 
later. This limiting form may be obtained by expressing Equation (2-45) in terms of its real 
and imaginary parts as 


a{P s } + { b) = a 2 (a8 + H){P S ) , 


2-46 


a{P 3 l = -(a8 + H){P„} - (1 + P){a} , 


where 


{P-} ={P*} +So{P a } . 2-48 

The column vectors {P^} and {P 3 } are the "stiffness" and "damping!' pressures respectively. 

If one formally sets o = 0 in Equation (2-46) it decouples from Equation (2-47) and may be 
solved directly. Since the right-hand side of Equation (2-46) becomes 0, the stiffness 
pressures are the same as those that would be obtained by computing the steady-state 
pressures at a perturbed eccentricity, subtracting the unperturbed pressures and dividing by the 
eccentricity perturbation. This latter method is frequently used for computing steady-state 
stiffnesses in incompressible flow and has been implemented here for the computation of 
"stiffnesses at 0 frequency" used in the above mentioned homing procedure. The 0 frequency 
damping pressures may be obtained by solving Equation (2-47) with {P^} as determined from 
the solution to Equation (2-46). 

The above discussion assumed that the perturbation coefficients in Equations (2-40) - (2-42) 
were determined prior to setting up the finite discretized equations for their solution. Identical 
results can be achieved by direct numerical perturbation of the difference equations. This 
approach, which has been implemented here and is described below, avoids algebraic error in 
determining the perturbation coefficients and may be used in complex situations where 
analytical determination of the perturbation coefficients is not feasible. 

After desired convergence of the Newton-Raphson process has been achieved under steady 
(unperturbed) conditions one may denote the resulting steady-state pressure vectors as { Pj } 
and the coefficient matrices as [C j ], etc. and Equation (2-34) may be written as 

[C'HPj} + [E']{P m } + [D']{P h } = {R 1 } . 2-49 

One may now perturb the kth component of the eccentricity vector by an amount rj, 
recalculate [C J ] at the new film thickness (but old pressure distribution, P) then subtract [C J ] at 
the old film thickness and divide the difference by r| to numerically obtain the derivative of 
[C j ] with respect to e k which will be denoted by [C j,k ]. Thus 
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The matrices [E jk ], [D i,k ] and { R j,k } are obtained in a similar' manner from the other coefficient 
matrices. If we introduce a disturbance to e k of magnitude e'r|, as was done in deriving 
Equation (2-43), then the change in the coefficient matrix [C j ] would be e'q[C j ’ k ] with 
corresponding changes in the other coefficient matrices. If we disturb Equation (2-49) by 
replacing {Pj} with {Pj} + rj{P' k }, [Cj] with [C J ] + e'q[C j,k ], etc. and collect terms of order 
r|, the following expression is obtained:If we set e' to unity in Equation (2-50) then {P' k } 


[C'HP'*} +[E i ]{P' | k _,} + [D i ]{P' J k „} =({ R i,k } -[C ik ]{P,} -[E |k ]{P M } -[D |k ]{P |tl })e' . 

2-50 


will become the 0 frequency stiffness pressure (the change in steadyrState pressure per unit 
change in eccentricity). It should be noted that the coefficients of the 0 frequency stiffness 
pressures in Equation (2-50) are the same as those for the steady-state pressures in Equation 
(2-49); only the right-hand side has changed. Equation (2-50) thus represents the construction 
of the discretized form of Equation (2-43) when o = 0. In order to complete the process for 
G * 0, one may introduce the same disturbances to the right-hand side of Equation (2-31), 
with H r = H + e'rj{a} and add the terms of order T| to the right-hand side of Equation (2-50). 
The terms to be added are -0([C j ]{Pj k }+{R j ’ k }e , )/5f where [C J ] are diagonal matrices whose 
components are 


Cj =(a8 + H) j+lj+1 AA 1 +(a8 +H) i+lj _ 1 AA 4 + (a8 +H) i _ Ij _ 1 AA 3 + (a8 +H) i _ lj+1 AA 2 


2-51 


and {R j,k } are column vectors whose components are 


R, |k =(1 ♦ p, AA, AA 4+ a i ;, 1 .,AA 3+ a i ;, h ,AA 2 )=(1 +P,)a,.]AA 


2-52 


The far right side of Equation (2-52) is a quadradically equivalent representation that was 
used in the computer program described in Section 3. One may now set e # = e 3at in 
Equation (2-50) and look for solutions in the form {Pj k } = {P* k }e 3m , by introducing these 
substitutions into Equation (2-50) and combining terms to obtain the final set of linear 
difference equations for the complex stiffness pressures { Pj k } : 


[C -iHPi*} +[E i ]{P j : i k } +[b i ]{P i : i k } ={R i,k } — [ O i,k ] { F> j } -[E ,,k ]{P H } -[D^HPjJ 


2-53 


where [C* J ] = [C J ] + 3 g[C j ] and {R J ’ k } = {R jk } - 3a{R J ’ k }. 
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The system of equations given by Equation (2-53) has been solved by the column method in a 
directly analogous manner to that used in solving Equation (2-34). The principal difference 
lies in the fact that all of the matrix operations were performed using complex arithmetic. 

The dimensionless, frequency dependent stiffness and damping coefficients were computed 
from the complex stiffness pressures in the previously described manner. Relationships of the 
following type may be used to calculate the physical stiffness and damping coefficients from 
the dimensionless ones: 


and 


K = K K 

rv xx rv O rv xx > 


K w =K 0 R 0 2 K 4 


k** = k 0 r 0 k x . 


2-54 


B xx = B o B xx 



=b 0 r 0 2 b„ 


B x* B 0 Ro B * 


2-55 


where Kq = p 0 R 0 2 /C and B 0 = 12 pRq 4 /C 3 . 


2.1.6 Optimization of Groove Parameters for Maximum Stagnation Pressure in a 
Concentric Cylindrical Seal 

Since spiral grooves are solely responsible for the axial stiffness of an aligned, gas-lubricated 
face seal with parallel surfaces under steady-state conditions, it is often desirable to optimize 
groove parameters for maximum axial stiffness. An optimization procedure for doing this has 
been implemented in the computer code SPIRALP described in Reference 1 1 . The analogous 
situation is not as evident in a concentric gas -lubricated cylindrical seal which will have 
considerable, if not maximum stiffness without spiral grooves. A large portion of the stiffness 
in the absence of spiral grooves will be cross coupled, particularly at low values of A, thus 
giving rise to stability problems which may be alleviated with the use of spiral grooves. The 
criteria for optimizing groove geometry from a dynamic standpoint would thus depend on 
both the desired load capacity and the various other elements in the system affecting 
rotordynamic performance. 

An alternate approach for developing a stand-alone criterion for optimizing groove geometry 
in a cylindrical seal is to maximize the pressure gradient that the grooves can generate at 
stagnation. If the grooves are being used to pump against a pressure gradient, the maximum 
stagnation pressure gradient would represent the maximum pressure gradient that the grooves 
could pump against without allowing any net flow to go through. It would also represent the 
maximum axial pressure gradient that the grooves could generate in an aligned, symmetric 
herringbone bearing in the absence of an imposed pressure gradient. In any event, the 
stagnation pressure gradient is a strong measure of spiral groove performance and even 
though optimizing it is not a precise criterion for optimizing dynamic performance, 
computations obtained with geometries optimized in this manner should provide a strong 
indication of the maximum benefits obtainable with the use of spiral grooves. 
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The stagnation pressure gradient for a cylindrical seal under concentric conditions may be 
obtained from Equation (2-22) by setting Q s = 0 (stagnation), 3P/30 = 0, H r = 1 (concentric), 
S = Z and R = 1 (cylindrical seal). The resulting equation may be solved for AP/AZ making 
use of the definition of A§ given by Equation (2-24) and the definitions of kj and k 4 given by 
Equation (2-25) to obtain the following relationship 

dP = A - 5ce(1 -(x)sinpcosp(r 3 -1) 
dZ “ a (l +I G ‘ 

The right-hand side of the above equation may be treated as a function of a, (3 and 
8(r =1+8) and has a maximum value of 3P/3Z = 0.09118 Ao) at a opt = 0.5, 

(3 opt = 0.2736 (15.68°) and 8 = 2.653. 

The variation of the pressure gradient near the optimum point is shown in Figure 4. The 
curve marked a was obtained by holding (3 and 8 at their optimum values and varying a. 

The other curves were obtained in an analogous manner. The curves show the sensitivity of 
the optimum pressure gradient to the various parameters and verify the existence of a relative 
maximum at the optimum point. 

Other approaches to the optimization problem are given in Reference 5 for spiral groove 
bearings and Reference 12 for spiral groove viscous pumps. 


2.2 Overview of Computer Code SPIRALG 

SPIRALG has been written to implement the analysis developed in the preceding section. 

Full input and output descriptions along with a number of sample runs are given in Reference 
13. A brief discussion of a few items not covered in Section 2.1 is given below. 

The analytical procedure contained in Section 2.1 has been oriented toward determining 
pressure distribution, load, flow, torque, stiffness and damping for a given film thickness 
distribution. In practice it is often desirable to determine the equilibrium film thickness or 
eccentricities from prescribed loads and possibly moments. SPIRALG provides a homing 
option for determining the eccentricities based on the steady-state bearing stiffnesses. This 
homing option is controlled by the input flag IHOME and is based on the procedure 
described below. 

If one were to write the dimensionless load and eccentricity as column vectors {W} and {£} 
(transpose row matrix [£]) and take the previous estimate (or initial guess) of {e} as {£} okl 
and the load vector computed from { e } old as (W} old , the steady-state stiffness matrix [K] 
could be used to arrive at a new approximation for {£}. The method for doing this is shown 
by first writing the equation for the change in load as {W} - (W} old = [K]({e} - (e} old ). The 
new approximation to {£} is obtained by inverting the stiffness matrix and solving for { e } as 
{£} = {£} old + [K] _ 1 ({W} - (W} old ). This approach is in effect the application of the Newton- 
Raphson method for determining the eccentricities. 
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While the above approach can be very effective it can also diverge if the initial guesses are 
bad. This divergence is usually accompanied by the generation of negative film thicknesses 
in the course of the iteration process. In order to attempt to correct this problem, an optional 
numerical damping algorithm has been implemented which replaces {£} with 
{£} = {e} old + p[K] _1 ({W} - {W} old ) when the originally calculated value of {£} would result 
in a negative film thickness. The numerical damping factor $ is determined by the program 
based on the value of the input parameter NUMDMP described in the input description. 

The cell method of discretization is designed to obtain quadratic accuracy. Numerical testing 
indicates that this has apparently been achieved. One may make use of this property to obtain 
greater accuracy, (or the same degree of accuracy with coarser grids and ensuing reductions in 
computer time) with the use of Romberg extrapolation. Suppose for example we computed 
the dimensionless torque T with a coarse grid and denoted it by T c then halved the grid 
spacing in both directions and recomputed T denoting it as T f (subscript denotes fine grid). If 
the truncation error were to approach 0 as the square of the grid spacing and T r were the true 
solution then (T f - T r ) = (T c - T r )/4, or T r = (4T f - T c )/3. The above extrapolation can, in 
principal, increase the rate of convergence from quadratic to cubic. The subroutine SPIRAL, 
provides the option of implementing Romberg extrapolation for all of the outputs by setting 
the flag IACC as explained in the input description. 

The logic used in SPIRALG for performing the pressure iterations, computing stiffness and 
damping coefficients, homing in on eccentricities and implementing Romberg extrapolation is 
shown in Figure 5. It can be seen there that when the homing process is implemented, it is 
completed for both coarse grid and fine grid solutions prior to performing the Romberg 
extrapolation. The extrapolation is thus performed with solutions obtained at two different 
displacements. When the displacements are specified (IHOME=0) extrapolations are 
performed with solutions obtained at the same displacement, which is believed to be a more 
accurate approach. If one were to compute displacements for a given loading and then 
recompute the loading from the displacements using Romberg extrapolation for both 
computations the computed loading would thus differ slightly from the input loading even 
though all tolerances were met. The degree of this difference will depend on the grid size 
and caution should be exercised in using Romberg extrapolation when homing in on the 
displacements with very coarse grids. 


2.3 Verification 

SPIRALG has been compared with the results of two computer codes. The first of these is 
MTI Computer Code PN471 which is fully described in Reference 6. The program is based 
on perturbation analyses and is only applicable to bearings and seals operating in the 
concentric, aligned position (£ x =£ y =(j>=\|/=0). The program does not predict loads and 
moments that would occur at finite displacements but it does predict stiffness and damping 
values as well as flow, torque, and power loss for spiral groove bearings as well as cylindrical 
and face seals. Since the program solves ordinary rather than partial differential equations it 
can be made to rapidly produce highly accurate results for evaluating the accuracy of 
SPIRALG. The second MTI computer code, named GASBEAR, is used to verify SPIRALG 
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Figure 5. Flow Diagram for Logic Used in SUBROUTINE SPIRAL 
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under displaced and misaligned conditions. GASBEAR was written for use in conjunction 
with plane journal bearings and cylindrical seals. It does not treat spiral grooves or face 
seals. Since SPIRALG does not contain and special instructions for treating concentric 
behavior and has relatively few instructions for distinguishing between face and cylindrical 
seals, the above two programs should provide reasonable verification. Since the treatment of 
the effects of spiral grooves under eccentric and misaligned conditions is believed to be new, 
terms that become significant only under those conditions remain unverified. 

The results of 7 verification tests are reported on the following pages. A SPIRALG output 
listing followed by the relevant output from the verification code, converted to equivalent 
units and format, is given for each case. The somewhat strange looking input values (unit 
ambient pressure, high RPM but low viscosity etc.) were selected to simplify the conversion 
process between dimensional quantities and the dimensionless ones that were used throughout 
the development of the code. The compressibility number of A=10 used for all cases and the 
seal pressure ratio of 2 used for imposed pressure gradients should be typical of many 
practical applications under fairly compressible conditions. 

Cases 1-6 show comparisons between SPIRALG and PN471. Romberg extrapolation was 
used for each of these cases, with 21 grid points in the circumferential direction and 5 sub- 
intervals in each of the two regions in the transverse direction for the coarse grid solution. 

The fine grid solutions thus use 41 circumferential points and 10 sub-intervals per region in 
the transverse direction. 

The first case verifies stiffness and damping values for a synchronous disturbance acting on a 
cylindrical seal with a herringbone groove pattern and an imposed pressure gradient. Cases 2 
and 3 verify the differences in stiffness and damping values predicted to occur for a 
cylindrical seal when grooves are placed on the rotor (with groove angles reversed) rather 
than on the stator. The static quantities (flow, torque and power loss) remain unchanged. 
Cases 4-6 show comparisons between SPIRALG and PN471 for a spiral groove face seal 
with an imposed pressure gradient at three different disturbance frequencies; zero (case 4), 
synchronous (case 5) and ten times synchronous (case 6). 

Case 7 shows a comparison between SPIRALG and GASBEAR for an eccentric, tilted 
cylindrical seal with an imposed pressure gradient. The grid size was chosen to match the 
maximum size allowable for the available version of GASBEAR. A separate program was 
written to perform Romberg extrapolations with the results of GASBEAR. The agreement 
between the two programs is good. The apparent discrepancy between the moments about the 
y axis is a result of the fact that the component is very small. The relative error obtained by 
dividing the discrepancy by the absolute magnitude of the moment vector is 0.33%. 
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( CASE 1 ) Concentric cyl . seal with pres. grad. A=10, a = 20 
SPIRAL GROOVE SHAFT SEAL, ROTATING SURFACE IS SMOOTH 

LENGTH, DIAMETER, CLEARANCE = 2.0000E+00, 2.0000E+00, 1.0000E-03 (IN) 

ROTATION SPEED, DISTURBANCE SPEED = 1.9099E+05, 1.9099E+05 (RPM) 

PRESSURE AT START, END AXIAL BOUNDARIES = 2.0000E+00, 1.0000E+00 (PSI) 

VISCOSITY - 8 . 3333E-11 (PSI-SEC), AMBIENT PRESSURE = 1.0000E+00 (PSI) 

CALCULATED FORCES IN X,Y DIRECTIONS = -1.0902E-15, -8.1153E-16 (LB) 
CALCULATED MOMENTS ABOUT X,Y AXES - -9.3565E-16, -2.4405E-16 (IN-LB) 

MINIMUM FILM THICKNESS = 1.0000E-03 (IN) 

FLOW = 1.2816E+01 (IN 3 /SEC) MEASURED AT 1.0000E+00 k (PSI) 

TORQUE = 1 . 7 015E- 02 (IN-LB), FILM POWER LOSS = 5.1562E-02 (HP) 

COMPRESSIBILITY NUMBER = 1.0000E+01, SQUEEZE NUMBER = .2.0000E+01 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP . UNIT ) 


DISP. 

x (IN) 

y (IN) 

0 (RAD) 

4r (RAD) 

FORCE UNIT 

Kx 

5 . 4607E+03 

9. 0379E+01 

-6 . 8886E+01 

2 . 6847E+02 

LB 

Ky 

-9 . 0379E+01 

5 . 4607E+03 

-2 . 6847E+02 

-6 . 8886E+01 

LB 

K0 

-1 . 0824E+02 

-1. 7773E+01 

7 . 0561E+02 

1 . 9521E+02 

IN- LB 

K0 

1 . 7773E+01 

-1 . 0824E+02 

-1 . 9521E + 02 

7 . 0561E+02 

IN-LB 

Bx 

7 . 3200E-02 

- 5 . 5832E- 02 

-1 . 5172E- 02 

-2 . 9538E-02 

LB-SEC 

By 

5 . 5832E- 02 

7 . 3200E-02 

2 . 9538E-02 

-1 . 5172E- 02 

LB-SEC 

B0 

1 . 9442E- 03 

6 . 9321E-03 

2 . 7091E- 02 

-1 . 1690E-02 

IN-LB-SEC 

Bl|/ 

-6 . 9321E- 03 

1 . 9442E- 0 3 

1 . 1690E- 02 

2 . 7091E-02 

IN-LB-SEC 



COMPARISON OF 

CASE 1 WITH 

PN4 71 


FLOW = 

1 . 282E+01 

(IN 3 /SEC) MEASURED AT 1.0000E+00 (PSI) 

TORQUE 

= 1.701E-02 

(IN-LB) , 

FILM POWER LOSS = 5.156E 

-02 (HP) 


DYNAMIC COEFFICIENTS ( 

FORCE UNIT / 

DISP. UNIT ) 


DISP. 

x (IN) 

y (IN) 

0 (RAD) 

0 (RAD) 

FORCE UNIT 

Kx 

5 . 4608E + 03 

9 . 0292E+01 

-6 . 8846E+01 

2 . 6823E + 02 

LB 

Ky 

-9 . 0292E+01 

5 . 4608E+03 

-2 . 6823E+02 

-6 . 8846E+01 

LB 

K0 

-1 . 0809E+02 

-1 . 7525E+01 

7 . 03 56E+02 

1 . 9537E+02 

IN- LB 

K0 

1 . 7525E + 01 

-1 . 0809E+02 

-1 . 9537E+02 

7 . 0356E+02 

IN-LB 

Bx 

7 . 3215E-02 

-5 . 5805E- 02 

-1 . 5192E-02 

-2 . 9534E- 02 

LB-SEC 

By 

5 . 5805E-02 

7 . 3215E-02 

2 . 9534E- 02 

-1 . 5192E- 02 

LB-SEC 

B0 

1 . 9403E-03 

6 . 9195E- 03 

2 . 7097E-02 

-1. 1667E-02 

IN-LB-SEC 

B0 

-6 . 9195E- 03 

1 . 9403E-03 

1 . 1667E- 02 

2 . 7097E-02 

IN-LB-SEC 
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( CASE 2 ) Concentric asymmetric cyl . seal, A=10, a^O 


SPIRAL GROOVE SHAFT SEAL, ROTATING SURFACE IS SMOOTH 

LENGTH, DIAMETER, CLEARANCE = 2.0000E+00, 2.0000E+00, 1.0000E-03 (IN) 

ROTATION SPEED, DISTURBANCE SPEED = 1.9099E+05, O.OOOOE+OO (RPM) 

PRESSURE AT START, END AXIAL BOUNDARIES = 1.0000E+00, 1.0000E+00 (PSI) 

VISCOSITY = 8 . 3333E-11 (PSI-SEC) , AMBIENT PRESSURE = 1.0000E+00 (PS 

CALCULATED FORCES IN X,Y DIRECTIONS = -2.8151E-15, -1.6233E-15 (LB) 
CALCULATED MOMENTS ABOUT X,Y AXES = 2.7524E-16, -5.2158E-17 ( IN- LB) 

MINIMUM FILM THICKNESS = 1.0000E-03 (IN) 

FLOW = 3 . 5000E + 00 (IN 3 /SEC) MEASURED AT 1.0000E + 00 (PSI) 

TORQUE = 1 . 8850E-02 (IN-LB), FILM POWER LOSS = 5.7122E-02 (HP) 

COMPRESSIBILITY NUMBER - 1.0000E+01, SQUEEZE NUMBER = . O.OOOOE+OO 


DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP . 

x (IN) 

y (IN) 

4> (RAD) 

(RAD) 

FORCE UNIT 

Kx 

4 . 96 57E + 03 

1 . 5689E + 03 

-4 . 9960E+02 

2 . 6673E + 02 

LB 

Ky 

-1 . 5689E+03 

4 . 9657E+03 

-2 . 6673E+02 

-4 . 9960E + 02 

LB 

K$ 

5 . 3012E+01 

-3 . 6877E+02 

6 . 0778E+02 

5 . 3984E+02 

IN-LB 

K4r 

3 . 6877E+02 

5 . 3012E + 01 

-5 . 3984E+02 

6 . 0778E+02 

IN-LB 

Bx 

-6 . 9953E-02 

- 1 . 1123E- 01 

4 . 3933E-02 

-1 . 2462E- 02 

LB-SEC 

By 

1 . 1123E- 01 

-6 . 9953E-02 

1 . 2462E- 02 

4 . 3933E-02 

LB-SEC 

B4> 

-8 . 3037E-03 

-4 . 4522E-03 

2 . 2471E- 02 

-4 . 3624E- 02 

IN-LB-SEC 

Bi)j 

4 . 4522E-03 

- 8 . 3037E- 03 

4 . 3624E- 02 

2 . 2471E-02 

IN-LB-SEC 



COMPARISON OF 

CASE 2 WITH 

PN4 71 


FLOW = 

3 . 500E+00 

(IN 3 /SEC) MEASURED AT 1.0000E+00 (PSI) 

TORQUE 

= 1 . 885E-02 

(IN-LB) , 

FILM POWER LOSS = 5.712E- 

02 (HP) 


DYNAMIC COEFFICIENTS ( 

FORCE UNIT / 

DISP. UNIT ) 


DISP. 

x (IN) 

y (IN) 

0 (RAD) 

(RAD) 

FORCE UNIT 

Kx 

4 . 9648E+03 

1 . 5692E+03 

-4 . 9982E+02 

2 . 6552E+02 

LB 

Ky 

-1 . 56 92E + 03 

4 . 9648E+03 

-2 . 6552E+02 

-4 . 9982E+02 

LB 

K<t> 

5 . 3156E + 01 

-3.6 793E + 02 

6 . 0602E+02 

5 . 3900E+02 

IN-LB 

KU 

3 . 6793E + 02 

5 . 3156E+01 

-5 . 3900E+02 

6 . 06 02E + 02 

IN-LB 

Bx 

-6 . 9955E-02 

-1 . 1143E- 01 

4 . 3975E-02 

-1 . 2344E-02 

LB-SEC 

By 

1 . 1143E-01 

-6 . 9955E- 02 

1 . 2344E- 02 

4 . 3975E-02 

LB-SEC 

Bcj> 

-8 . 3373E-03 

-4 . 4940E- 03 

2 . 2617E-02 

-4 . 3536E-02 

IN-LB-SEC 

BV 

4 . 4940E- 03 

-8 . 3373E-03 

4 . 3536E-02 

2 . 2617E-02 

IN-LB-SEC 
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( CASE 3 ) Same case with grooves on moving surf., groove angle rev. 


SPIRAL GROOVE SHAFT SEAL, ROTATING SURFACE IS GROOVED 

LENGTH, DIAMETER, CLEARANCE = 2.0000E+00, 2.0000E+00, 1.0000E-03 (IN) 

ROTATION SPEED, DISTURBANCE SPEED = 1.9099E+05, O.OOOOE+OO (RPM) 

PRESSURE AT START, END AXIAL BOUNDARIES = 1.0000E+00, 1.0000E+00 (PSI) 

VISCOSITY - 8.3333E-11 (PSI-SEC), AMBIENT PRESSURE = 1.0000E+00 (PSI) 

CALCULATED FORCES IN X,Y DIRECTIONS = -3.4729E-16, -1.2658E-15 (LB) 

CALCULATED MOMENTS ABOUT X,Y AXES = 1.6213E-16, 6.2386E-17 (IN-LB) 

MINIMUM FILM THICKNESS = 1.0000E-03 (IN) 

FLOW = 3 . 5000E+00 (IN 3 /SEC) MEASURED AT 1.0000E+00 ‘(PSI) 

TORQUE = 1 . 8850E- 02 (IN-LB), FILM POWER LOSS = 5.7122E-02 (HP) 

COMPRESSIBILITY NUMBER = 1.0000E+01, SQUEEZE NUMBER = „ 0.0000E+00 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP . 

x (IN) 

y (IN) 

0 (RAD) 

l]/ (RAD) 

FORCE UNIT 

Kx 

4 . 5241E+03 

1 . 2 933E+03 

-3 . 6977E + 02 

2 . 1989E+02 

LB 

Ky 

-1 . 2933E+03 

4 . 5241E+03 

-2 . 1989E+02 

-3 . 6977E+02 

LB 

K0 

3 . 4694E+02 

-4 . 7068E+02 

5 . 9016E+02 

5 . 1362E+02 

IN-LB 

Kijr 

4 . 7068E+02 

3 . 46 94E+02 

-5 . 1362E+02 

5 . 9016E+02 

IN-LB 

Bx 

-2 . 3432E-02 

-1.2667E-01 

3 . 4823E-02 

-5 . 9373E- 03 

LB-SEC 

By 

1 . 2667E-01 

-2 . 3432E- 02 

5 . 9373E-03 

3 . 4823E- 02 

LB-SEC 

Bcj) 

-2 . 7769E- 02 

6 . 0874E-03 

2 . 3 559E-02 

-4 . 4038E-02 

IN-LB-SEC 

Blp- 

-6 . 0874E-03 

-2 . 7769E- 02 

4 . 4038E-02 

2 . 3559E-02 

IN-LB-SEC 



COMPARISON OF 

CASE 3 WITH 

PN4 71 


FLOW = 

3 . 500E+00 

(IN 3 /SEC) MEASURED AT 1.0000E+00 (PSI 

) 

TORQUE 

= 1.885E-02 

(IN-LB) , 

FILM POWER LOSS = 5.712E- 

02 (HP) 


DYNAMIC COEFFICIENTS ( 

FORCE UNIT / 

DISP. UNIT ) 


DISP . 

X (IN) 

y (IN) 

0 (RAD) 

)\! (RAD) 

FORCE UNIT 

Kx 

4 . 5236E+03 

1 . 2934E+03 

-3 . 6992E+02 

2 . 1883E+02 

LB 

Ky 

-1 .2934E + 03 

4 . 5236E+03 

-2 . 1883E+02 

-3 . 6992E+02 

LB 

Kef) 

3 . 46 97E+02 

-4 . 6939E+02 

5 . 8 843E + 02 

5 . 12 92E + 02 

IN-LB 

Kij ; 

4 . 6939E+02 

3 . 4697E+02 

-5 . 1292E+02 

5 . 8843E+02 

IN-LB 

Bx 

-2 . 3 4 7 IE - 02 

-1 . 2673E-01 

3 . 4849E- 02 

-5 . 8520E- 03 

LB-SEC 

By 

1 . 2673E-01 

-2 . 3471E- 02 

5 . 8 52 0E- 03 

3 . 4849E- 02 

LB-SEC 

B4> 

-2 . 7798E-02 

5 . 9155E-03 

2 . 3693E-02 

-4 . 3962E-02 

IN-LB-SEC 

B\J/ 

- 5 . 9155E-03 

-2 . 7798E-02 

4 . 3962E-02 

2 . 3693E-02 

IN-LB-SEC 
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( CASE 4 ) Face seal, no tilt, with pres, grad., A=10, a=0 
SPIRAL GROOVE FACE SEAL, ROTATING SURFACE IS SMOOTH 

ID, OD, REFERENCE FILM THICKNESS = 1.0000E+00, 2.0000E+00, 1.0000E-03 (IN) 

ROTATION SPEED, DISTURBANCE SPEED = 1.9099E+05, O.OOOOE+OO (RPM) 

INSIDE, OUTSIDE PRESSURE = 2.0000E+00, 1.0000E+00 (PSI) 

VISCOSITY = 8 . 3333E-11 (PSI-SEC) , AMBIENT PRESSURE = 1.0000E+00 (PSI) 

CALCULATED FORCE IN Z DIRECTION = 1.3612E+00 (LB) 

CALCULATED MOMENTS ABOUT X,Y AXES = -3.0782E-16, 1.1794E-15 ( IN-LB) 

MINIMUM FILM THICKNESS = 1.0000E-03 (IN) 

FLOW = 2 . 2763E+01 (IN 3 /SEC) MEASURED AT 1.0000E+00 ‘ (PSI) 

TORQUE = 2 . 2645E- 03 (IN-LB), FILM POWER LOSS = 6.8621E-03 (HP) 

COMPRESSIBILITY NUMBER = 1.0000E+01, SQUEEZE NUMBER = .0.0000E+00 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP. 

Z (IN) 

cj> (RAD) 

i|r (RAD) 

FORCE UNIT 

Kz 

3 . 9524E+02 

-4 . 8585E-04 

-4 . 8576E-04 

LB 

K(J> 

-6 . 1349E-07 

2 . 0998E + 02 

8 . 2256E + 01 

IN-LB 

Kty 

-3 . 0546E-07 

-8 . 2256E+01 

2 . 0998E+02 

IN -LB 

Bz 

2 . 9491E-02 

-7 . 3530E-07 

- 7 . 3551E-07 

LB-SEC 

B0 

-1 . 12 8 0E- 09 

7 . 0914E- 03 

-1 . 8082E-03 

IN-LB-SEC 

Blj; 

-2 . 3148E- 0 9 

1 . 8082E- 03 

7 . 0914E-03 

IN-LB-SEC 


COMPARISON OF CASE 4 WITH PN471 


CALCULATED FORCE IN Z DIRECTION = 1.3612E+00 (LB) 


FLOW = 

2 . 2763E+01 

(IN 3 /SEC) MEASURED AT 1 

. 0000E+00 (PSI) 

TORQUE 

= 2 . 2644E- 03 

(IN-LB) , 

FILM POWER 

LOSS = 6 . 8620E-03 


DYNAMIC COEFFICIENTS ( 

FORCE UNIT / 

DISP. UNIT ) 

DISP. 

z (IN) 

cj) (RAD) 

4r (RAD) 

FORCE UNIT 

Kz 

3 . 9524E+02 



LB 

K4> 


2 . 0998E+02 

8 . 2240E+01 

IN-LB 

K\|r 


-8 . 2240E+01 

2 . 0998E+02 

IN- LB 

Bz 

2 . 9485E-02 



LB-SEC 

B4> 


7 . 0901E-03 

-1. 8058E-03 

IN-LB-SEC 

Bij; 


1 . 8058E-03 

7 . 0901E- 03 

IN-LB-SEC 
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( CASE 5 ) Face seal, no tilt, with pres, grad., A=10, o=20 


SPIRAL GROOVE FACE SEAL, ROTATING SURFACE IS SMOOTH 

ID, OD, REFERENCE FILM THICKNESS = 1.0000E+00, 2.0000E+00, 1.0000E-03 (IN) 

ROTATION SPEED, DISTURBANCE SPEED = 1.9099E+05, 1.9099E+05 (RPM) 

INSIDE, OUTSIDE PRESSURE = 2.0000E+00, 1.0000E+00 (PSI) 

VISCOSITY = 8 . 3333E-11 (PSI-SEC) , AMBIENT PRESSURE = 1.0000E+00 (PSI) 

CALCULATED FORCE IN Z DIRECTION = 1.3612E+00 (LB) 

CALCULATED MOMENTS ABOUT X,Y AXES = -3.0782E-16, 1.1794E-15 (IN-LB) 

MINIMUM FILM THICKNESS = 1.0000E-03 (IN) 

FLOW = 2 . 2763E+01 (IN 3 /SEC) MEASURED AT 1.0000E+00 ‘(PSI) 

TORQUE = 2 . 2645E- 03 (IN-LB), FILM POWER LOSS = 6.8621E-03 (HP) 

COMPRESSIBILITY NUMBER = 1.0000E+01, SQUEEZE NUMBER = ,2.0000E+01 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP. 

z (IN) 

4> (RAD) 

ij; (RAD) 

FORCE UNIT 

Kz 

5 . 3566E+02 

-2 . 9756E- 04 - 

2 . 9746E- 04 

LB 

Kp 

-2 . 1294E- 07 

2 . 4013E+02 

7 . 0968E+01 

IN -LB 

K\J/ 

2 . 5103E-07 

-7 . 0968E+01 

2 . 4013E+02 

IN-LB 

Bz 

2 . 7697E-02 

3 . 8068E-09 

3 . 8070E-09 

LB-SEC 

B0 

2 . 2173E-12 

6 . 7533E-03 - 

1 . 6346E- 03 

IN-LB-SEC 

BiJ/ 

-2 . 0783E-12 

1 . 6346E-03 

6 . 7533E- 03 

IN-LB-SEC 



COMPARISON OF 

CASE 5 WITH 

PN471 


CALCULATED FORCE IN Z DIRECTION = 1.3612E+00 (LB) 

FLOW = 2 . 2763E + 01 (IN 3 /SEC) MEASURED AT 1.0000E+00 (PSI) 

TORQUE = 2 . 2644E-03 (IN-LB), FILM POWER LOSS = 6.8620E-03 (HP) 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 

FORCE UNIT 
LB 

IN- LB 
IN-LB 
LB-SEC 
IN-LB-SEC 
IN-LB-SEC 


NAS A/CR— 2004-2 1 3 1 99/VOL3 34 


DISP. 

z (IN) 

()) (RAD) 

4/ (RAD) 

Kz 

5 . 3557E + 02 



K<(> 


2 . 4 0 1 1E + 02 

7 . 0959E + 01 

KtJ; 


-7 . 0959E+01 

2 . 4011E + 02 

Bz 

2 . 7692E-02 



B4> 


6 . 7528E-03 

-1 .6341E-03 

Bd 


1 . 6341E-03 

6 . 7528E-03 



( CASE 6 ) Face seal, no tilt, with pres, grad., A=10, o - 200 
SPIRAL GROOVE FACE SEAL, ROTATING SURFACE IS SMOOTH 

ID, OD, REFERENCE FILM THICKNESS = 1.0000E+00, 2.0000E+00, 1.0000E-03 (IN) 

ROTATION SPEED, DISTURBANCE SPEED = 1.9099E+05, 1.9099E+06 (RPM) 

INSIDE, OUTSIDE PRESSURE = 2 . 0000E+00 , 1.0000E+00 (PSI) 

VISCOSITY = 8 . 3333E-11 (PSI-SEC), AMBIENT PRESSURE = 1.0000E+00 (PSI) 

CALCULATED FORCE IN Z DIRECTION - 1.3612E+00 (LB) 

CALCULATED MOMENTS ABOUT X,Y AXES = -3.0782E-16, 1.1794E-15 (IN-LB) 

MINIMUM FILM THICKNESS = 1.0000E-03 (IN) 

FLOW = 2 . 2763E+01 (IN 3 /SEC) MEASURED AT 1.0000E+00 ‘(PSI) 

TORQUE = 2 . 2645E-03 (IN-LB), FILM POWER LOSS = 6.8621E-03 (HP) 

COMPRESSIBILITY NUMBER = 1.0000E+01, SQUEEZE NUMBER = . 2.0000E+02 


DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP . 

z (IN) 

4> (RAD) 

0 (RAD) 

FORCE UNIT 

Kz 

2 . 3872E+03 

-4 . 1732E-05 

-4 . 1685E-05 

LB 

Kp 

-2 . 5267E- 08 

7 . 2912E+02 

- 7 . 4833E+00 

IN- LB 

Kf 

5 . 8389E-08 

7 . 4833E+00 

7 . 2912E + 02 

IN-LB 

Bz 

4 . 1358E-03 

5 . 3850E-10 

5 . 3818E-10 

LB-SEC 

B0 

4 . 2377E-13 

1 . 2426E- 03 

-5 . 8922E-05 

IN-LB-SEC 

Bi]j 

-5 . 5928E-13 

5 . 8922E-05 

1 . 2426E-03 

IN-LB-SEC 


COMPARISON OF CASE 6 WITH PN471 
CALCULATED FORCE IN Z DIRECTION = 1.3612E+00 (LB) 

FLOW = 2.276 3 E + 0 1 (IN 3 /SEC) MEASURED AT 1.0000E+00 (PSI) 

TORQUE - 2 . 2 6 44E- 03 (IN-LB), FILM POWER LOSS = 6.8620E-03 (HP) 


DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP . 

z (IN) 

0 (RAD) 

0 (RAD) 

FORCE UNIT 

Kz 

2 . 3870E + 03 



LB 

K0 


7 . 2907E+02 

-7 . 4987E + 00 

IN- LB 

K0 


7 . 4987E+00 

7 . 2907E+02 

IN-LB 

Bz 

4 . 1383E-03 



LB-SEC 

B0 


1 . 2436E-03 

-5 . 8939E-05 

IN-LB-SEC 

Bi(i 


5 . 8939E-05 

1 . 2438E- 03 

IN-LB-SEC 
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( CASE 7 ) Misaligned shaft seal with pres. grad, to comp, w/ GASBEAR 
SPIRAL GROOVE SHAFT SEAL, ROTATING SURFACE IS SMOOTH 

LENGTH, DIAMETER, CLEARANCE = 2.0000E+00, 2.0000E+00, 1.0000E-03 (IN) 

ROTATION SPEED, DISTURBANCE SPEED = 1.9099E+05, O.OOOOE+OO (RPM) 

PRESSURE AT START, END AXIAL BOUNDARIES = 2.0000E+00, 1 . 0000E+00 (PSI) 

VISCOSITY = 8 . 3333E- 11 (PSI-SEC), AMBIENT PRESSURE = 1.0000E+00 (PSI) 

CALCULATED FORCES IN X,Y DIRECTIONS = 3.1535E+00, 1.2350E+00 (LB) 

CALCULATED MOMENTS ABOUT X,Y AXES - 3.7793E-01, -1.9883E-03 (IN-LB) 

MINIMUM FILM THICKNESS = 5.0000E-04 (IN) 

FLOW - 4 . 7315E+00 (IN 3 /SEC) MEASURED AT 1.0000E + 00 dPSI) 

TORQUE = 2 . 3248E- 02 (IN-LB), FILM POWER LOSS = 7.0448E-02 (HP) 


COMPRESSIBILITY NUMBER = 1.0000E+01, SQUEEZE NUMBER = . 0.0000E+00 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP. 

x (IN) 

y (IN) 

$ (RAD) 

l|/ (RAD) 

FORCE UNIT 

Kx 

9 . 2287E+03 

4 . 8183E+03 

1 . 1084E+03 

5 . 7870E + 02 

LB 

Ky 

-2 . 6167E+03 

9 . 5964E + 03 

4 . 7033E + 02 

8 . 7124E+02 

LB 

Kc)> 

- 9 . 8594E + 01 

6 . 5068E+02 

1. 5189E+03 

7 . 7760E+02 

IN-LB 

Klfr 

7 . 1674E+02 

-9 . 7210E+01 

-1 . 0491E+03 

1 . 2765E + 03 

IN-LB 

Bx 

-2 . 0553E- 01 

-4 . 5816E-01 

-1. 0900E-01 

6 . 2037E-02 

LB -SEC 

By 

2 . 5253E-01 

-2 . 9852E- 01 

-6 . 7612E- 02 

-7 . 7454E- 02 

LB- SEC 

B$ 

2 . 3869E-02 

1 . 4706E- 02 

3 . 7737E-02 

-8 . 3553E-02 

IN-LB-SEC 

B4r 

- 7 . 4893E- 02 

-5 . 8155E- 02 

5 . 7003E- 02 

6 . 7719E- 02 

IN-LB-SEC 


COMPARISON OF CASE 7 WITH GASBEAR USING SAME GRID AND ROMBERG EXTRAPOLATION 


CALCULATED FORCES IN X,Y DIRECTIONS = 3.1529E+00, 1.2362E+00 (LB) 

CALCULATED MOMENTS ABOUT X,Y AXES = 3.7751E-01, -7.3904E-04 (IN-LB) 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP . 

x (IN) 

y (IN) 

<t> (RAD) 

(RAD) 

FORCE UNIT 

Kx 

9 . 2287E+03 

4 . 8097E + 03 

1 . 1044E + 03 

5 . 7867E+02 

LB 

Ky 

-2 . 6113E+03 

9 . 5993E+03 

4 . 7249E+02 

8 . 7343E+02 

LB 

K$ 

-1 . 0249E+02 

6 . 5339E+02 

1 . 5237E+03 

7 . 7 076E + 02 

IN-LB 

Kijr 

7 . 14 55E + 02 

- 9 . 4976E+01 

-1 . 0403E+03 

1 . 2793E + 03 

IN- LB 

Bx 

-2 . 0 5 6 IE- 01 

-4 . 5802E-01 

-1 . 0909E-01 

6 . 1715E- 02 

LB -SEC 

By 

2 . 5292E-01 

-2 . 9885E- 01 

-6 . 7494E-02 

-7 . 7535E-02 

LB-SEC 

B0 

2 . 3875E-02 

1 . 4639E-02 

3 . 7502E-02 

-8 . 3706E-02 

IN-LB-SEC 

Bi]/ 

-7 . 4561E-02 

-5 . 7965E-02 

5 . 6982E-02 

6 . 7347E-02 

IN-LB-SEC 
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3.0 INDUSTRIAL CODE SPIRALI - INCOMPRESSIBLE, 
TURBULENT SPIRAL-GROOVED CYLINDRICAL 
FACE SEALS 


Spiral groove bearings and seals are used to provide stability, load support and pumping for both 
cylindrical and face seal geometries. In the case of a cylindrical bearing, grooves are usually 
designed to pump against each other in a symmetric arrangement to provide enhanced stability. A 
lightly loaded cylindrical seal operating at a low axial flow rate will produce a force that is nearly 90 
degrees out of phase with the displacement which will tend to destabilize the rotating shaft. The 
introduction of spiral grooves can significantly increase the component of force in phase with the 
displacement and decrease the out of phase component thereby improving stability. 

In the case of a face seal or thrust bearing, spiral grooves are often introduced as the primary 
means of load support. Since a symmetric arrangement is not possible in a radial geometry, the 
grooves are usually designed to pump towards an ungrooved dam region. The resistance of the 
dam region increases as the film thickness decreases hence the pumping pressure rise increases 
thereby giving rise to a positive axial stiffness. The spiral grooves can also be used to pump 
against an applied pressure gradient thereby resulting in either reduced or reversed leakage. 

The computer code SPIRALI predicts performance characteristics of incompressible cylindrical and 
face seals. Performance characteristics include load capacity (for face seals), leakage flow, power 
requirements and dynamic characteristics in the form of stiffness, damping and apparent mass 
coefficients in 4 degrees of freedom for cylindrical seals and 3 degrees of freedom for face seals. 
These performance characteristics are computed as functions of seal and groove geometry, load or 
film thickness, running and disturbance speeds, fluid viscosity, and boundary pressures. 

The basic assumptions that have gone into the computer code are listed below: 

1 . The flow is assumed to be isothermal and incompressible. 

2. Turbulence is treated with an extended form of the Hirs bulk flow model (Reference 14), 
generalized to include separate and arbitrary friction factor - Reynolds number 
relationships for each surface. 

3. Inertia effects (which throughout this report will refer to additional effects to those 
inherent in turbulence) associated with film discontinuities will be treated with the use of 
loss coefficients. 

4. Circumferential and transient effects are treated with the use of small perturbations to a 
steady state first-order solution for a concentric, aligned seal. These effects are 
characterized by either by stiffness and damping coefficients that are dependent on the 
disturbance frequencies or by stiffness, damping, and apparent mass coefficients. 

5. The film thickness is assumed to be small compared with seal lengths and diameters but 
large compared with surface roughness. 
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6. Narrow groove theory is used to characterize the effects of grooves by a global pressure 
distribution without requiring computations on a groove by groove basis. This has 
previously involved neglecting edge effects and local inertia effects associated with 
groove to groove pressure variations. In general, narrow groove theory is valid when 
there are a sufficiently large number of grooves so that 2Tt/N g « 1, where N g is the 
number of grooves. This theory has been extended to account for the pressure jumps 
arising from the effects of local inertia associated with the sudden changes in cross 
section at the inlet to each groove and ridge which can be of the order of N g larger than 
the other local inertia effects which have been neglected. 

7. Although transverse (axial for cylindrical seals or radial for face seals) variations in the 
film thickness profdes are permitted, machined surfaces for seals are assumed to be 
axisymmetric with the exception of effects of spiral grooves. 

The above assumptions still leave the code applicable to a broad range of cylindrical and face 
seal applications. Practical designs should contain a fairly large number of grooves to ensure 
smooth, isotropic operation. Elastic and thermal distortions as well as machining tolerances 
should also be estimated to validate the axisymmetric clearance assumption and may be 
characterized in part by inclusion of film variations in the transverse direction. The overall 
accuracy of the program will depend on the grid size used. Factors such as large film thickness 
variations in the transverse direction, inclusion of very small transverse inertia effects and large 
values of the length to diameter ratio could require an increased number of grid points. 

A derivation of the equations governing the performance of turbulent, incompressible, spiral 
groove cylindrical and face seals along with a description of their solution is given in the next 
section. This will be followed by a description of the computer codes including an input 
description, sample cases and comparisons with results of other codes. 

3.1 Formulation and Method of Solution 

Coordinate variables will be used to make the equations developed here applicable to both 
cylindrical and face seals and may be defined with the aid of Figure 6. The circumferential 
coordinate, 0, is as shown in Figure 6. The transverse coordinate is described by the variable, s, 
which is taken to equal the radial coordinate, r, for a face seal and the axial coordinate, z, for a 
cylindrical seal. The quantity r, when it appears will denote radial position for a face seal and 

should be set equal to the shaft radius, R 0 for a cylindrical seal. The symbols i , j will be used 
to denote unit vectors in the 0 and s directions respectively. 


Formulation of Bulk Flow Equations for Turbulent Lubrication 

The equations for turbulent bulk flow to be presented here will be based on a generalization of 
the Hirs model. Figure 7 shows the velocities and forces associated with a differential element 
in the 0 direction. The bulk flow integrated momentum equations for the 0 and s directions are 


/ 

ph 

V 


0U du u gu 

dt 0s r 89 


uvl f 

r ; 


r 00 


3-1 
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Figure 6. Coordinate System for Spiral Groove Analysis 
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The integrated continuity equation is 


- — (rvh).i— (uh).— - 0 
r as r ae at 


3-3 


In the above equations u and v represent the components of the velocity vector, u=u i + vj,p 

denotes the pressure, x a and T b are the shear stresses acting on the fluid at surfaces (a) and (b) 
respectively as shown in Figure 7 and I f is a flag parameter set equal to 1 for a face seal and 0 for 
a cylindrical seal. 

Each of these shear stresses is assumed to act in the direction of the fluid velocity relative to the 
respective surface and is related to that velocity by means of a coefficient of resistance or 
friction factor as follows: 


\ I U - u a I f a 


2hp | u - u. 


fu-u a) - ■ 


for surface (a) and for surface (b) 


\ ' -^Plu A! f b 


/ 2hp|u-u b j > 


(U-U b ) ■ -VT R bfb( R b )(u-U b) . 
4 h 


3-4 


3-5 


where the shear factors f a and f b are one fourth of the resistance coefficients or friction factors 
for the respective surfaces (Reference 15) and are functions of the respective Reynolds numbers 
R a and R b which are in turn defined as 

R a ■ 2h | u - u a | p/p , R b - 2h | u - u b | p/p . 3-6 

The use of various shear functions in the prediction of turbulent behavior of bearings and seals 
may be found in References 16 - 19. Power law relationships of the form f = nR m have been 
used in References 16 - 18 to characterize both smooth and rough surfaces. If both surfaces are 
smooth the relationship 

f a (RJ " n o R a m ° . f b (R b ) - n 0 R b m ° 3-7 

may be used with n 0 and m 0 obtained from the Blasius Equation (Reference 14) as n 0 = 0.0751, 
mo = -.25 for turbulent flow (Reynolds numbers ranging from 4,000 to 100,000). The 
corresponding values for laminar flow are n 0 = 24 and m 0 = -1 . 
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Equations (3-1) - (3-3) may be expressed in dimensionless form with the use of the following 
dimensionless variables 


h- 


h 

C ' 



S- 



f 


^OT 

C Po 


and parameters 


2CV 0 P 

M 


R 


2C 


R 



. P- 


P V Q r D 

4C 2 p 0 ‘ 


3-8 


3-9 


In the above equations V 0 represents a characteristic velocity, p 0 a characteristic pressure, r 0 and 
C represent the radius and clearance for a cylindrical seal and the outside radius and 
characteristic film thickness for a face seal. 

Henceforth surface (b) will be taken to be stationary and surface (a) will be taken to move with a 
circumferential velocity ro). The Reynolds numbers R a and R b given by Equation (3-6) then become 


R a - Rh\/(u-ftS) 2 + v 2 , R b -Rhv/u 2 + v 2 . 3-10 

The integrated momentum and continuity equations may now be expressed in the dimensionless 
form 


1 dp 
r 30 


pr<D(u,v,h) 


/ 

prR* 

V 


3u - 3u 

at ds 


u 3G | | uv ^ 
r 30 f r 


3p 

as 


prMJ(u f v,h) 


prR* 


( - 
3v 

, at 



u av 

r as 



3-11 


3-12 


where 


and 


1 a(rvh) + 1 a(uh) 3 h 

f as f a0 


0(u,v,h) = 


(u-f(S)R a f a (RJ + uR b f b (R b ) 
h 2 


MJ(u,v,h) - 


R a f a(RaMVb(Rb> , 

h 2 


V 


3-13 


3-14 


3-15 
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Equations (3-1 1) - (3-13) together with the functions defined by Equations (3-10), (3-14) and (3- 
15) constitute the dimensionless differential equations for the generalized Hirs bulk flow model 

under consideration. These equations require a prescribed film thickness profile, h(S,0, t), 

definitions of f a and f b such as those given by Equation (3-7) and boundary conditions which will 
be treated below. 


Boundary and Continuity Conditions 

This analysis is limited to essentially axisymmetric surface shapes and constant boundary 
pressures. Although spiral grooves will be considered, the groove angle, depth and spacing will 
be independent of 0. Individual pads and partial films are not treated here. As such, all 0 (and 
time) dependence will be associated with displacements of the sealing surfaces. The only 
conditions on 0 are the requirement of periodicity. That is, pressure and flow rates at 0 + 2 tx 
must equal those at 0 for all 0. 

Variations of surface shapes in the s direction, including discontinuities are considered. The 
upstream values of v and p at discontinuities in h are denoted by v ; and pj as shown in Figure 8. 
The film thickness changes from the upstream value of h - Ah to the downstream value of h 
across the jump. 

Since the flow rate normal to the jump must be continuous 


(h-Ah)Vj - hv @ s - Sj , 
where Sj denotes the s coordinate at the jump. 

The pressure change at s ; may be expressed in terms of a loss coefficient based on 
downstream conditions (Reference 20) as follows: 

Pj * |p v j 2 - P * ~pv 2 (1 *5) @ S - Sj . 


p , v , h - Ah 

j j 

p, v, h 

-Ah 

s flow direction 


95TR34-V3 


Figure 8. Jump in Film Thickness 
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The preceding two equations may be rearranged and expressed in dimensionless form as 


V . (1 -f)ij , Ap - p - Pj - |prR-x(v,h,Ah) @S - Sj , 3 _ 16 

h ^ 

where 




/ - \ 

2 

X(v,h,Ah) - - 

1 - 

h 

+ ^ 



v h - Ah j 



The loss coefficient due to a contraction will in general be a function of the axial Reynolds 

number and will be denoted by the function ((Rhv) but is frequently taken to be constant 
(References 16 - 18). The loss coefficient for an expansion may be found in Reference 20 and is 
given together with the corresponding relationship for a contraction below: 


C(R,h,v) 

( 


1 


l V 


h-Ah 


, Ah < 0 (contraction) 
, Ah > 0 (expansion) 


3-18 


The above considerations together with the requirement that the tangential component of 
velocity is continuous at jumps across constant s boundaries provides sufficient information to 
prescribe the boundary conditions at the inlet and exit to the seal as well as the continuity 
conditions at jumps. 

The tangential component of the velocity, u in and the pressure, p in are prescribed at the inlet side 
of the seal. The relationship between p in and v at the inlet is obtained by setting p ; = p jn and 

evaluating % as Ah -> -<» in Equations (3-16) - (3-18). These conditions may be written in 
dimensionless form as 


a ■ Q in ■ P ■ Pin * . p. n _ iprFftl.Qv 2 @S-S ln . 3 _ 19 

The pressure at the exit from the seal will be taken to be continuous and equal to the exit side 
pressure, p ex . 


P - Pex @ s - s ex - 3-20 

It should be noted that S in may correspond to either the left or right hand sides of the seal (lower 
or upper values of S) depending on the direction of flow. S ex will be the value of S at the 
opposite side of the seal from S = S in . 
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Narrow Groove Theory for Spiral Groove Seals 

The development of the narrow groove equations to be presented here will expand upon 
equations used in Reference 21 which, in turn, contains a reformulation of the lubrication 
equations for laminar, compressible spiral groove seals presented in Reference 22. 

The purpose of narrow groove theory is to obtain equations and solutions for continuous, smooth 
global pressure and velocity distributions that approximate the limiting form of the solution to 
Equations (3-10) - (3-15) as the number of grooves becomes large, with the groove angle, p and 
the groove to pitch ratio, oc, held constant. This process greatly reduces the complexity of the 
solution and eliminates the need for the large grid that would be required to adequately describe 
all of the grooves individually. 

The discontinuities in film thickness associated with the grooves will give rise to discontinuities 
in the pressures, pressure gradients, and velocities at the ridge-groove interfaces. The relative 
discontinuities in the velocity components normal to the grooves will be of the same order as the 
relative film discontinuities. The predicted discontinuities in the local pressures at the interfaces 
would be of the order of the inertia (p*R*) terms on the right hand sides of Equations (3-11) and 
(3-12), It is important to note however, that as both the number of grooves and the inertia terms 
become large, flows will not be able to fully develop within or between grooves. 

Figure 9 shows three conceptual states of flow development in the order of increasing 
significance of local fluid inertia and decreasing degree of flow development. These are Case a, 
where the local pressure changes at the inlet to each groove (or ridge) are negligible; Case b, 
where a wake forms at the inlet to a groove whose extent is small compared to the length of a 
groove and Case c, where a jet forms and prevails over most or all of the length of the groove. 
These three cases are discussed individually below: 

a. Negligible inlet pressure changes. Spiral groove seals frequently operate at relatively 
small clearances (as opposed to labyrinth seals or even damping seals) and will pump at 
relatively low flow rates compared with those produced by imposed pressure gradients of 
sufficient magnitude to produce large inertia effects. Under these conditions, the flow 
can be treated as fully developed over the entire groove and inertia effects may be 
neglected in the analysis of the local pressures and velocities. All inertia effects in 
Equations (3-11) - (3-20) may, however, be included in the treatment of the global 
pressure and velocities. As such, pumping affected by turbulent shear in regions 
containing spiral grooves would be accounted for and global inertia effects due to 
circular steps and grooves and other axial clearance variations would still be included in 
all regions. 

A schematic showing a global and local pressure distribution for this case along with 
some of the spiral groove film variables under consideration is presented in Figure 10. 

The local pressure profile is shown by the sawtooth lines whose lower vertices, for 
purposes of illustration, are connected by the global pressure profile, p. The global 
pressure profile does not necessarily lie at the lower vertices of the local pressure profile 
but could lie anywhere between the lower and upper vertices. In the limit as the number 
of grooves becomes large, curves connecting the upper and lower vertices will approach 
each other (subject to the approximations discussed above). This limiting behavior is not 
true of u, v, h or the local pressure gradients which will have different values over lands 
and grooves no matter how large the number of grooves. 
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Figure 9. Examples of Degrees of Film Formation in Groove 
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b. Fully developed films with significant inlet pressure changes. This case extends Case 
a by including local pressure changes at the inlet to each groove and ridge with the use of 
loss coefficients in a manner analogous to that used at the inlets to each region given by 
Equations (3-16) - (3-18). The treatment of these pressure changes as boundary pressure 
discontinuities assume that the length of the wake shown in Figure 9b is small compared 
with other relevant circumferential and transverse length dimensions. The relevant 
length for comparison with the wake length at the start of a region would be the 
transverse length of the region. In the case of spiral grooves within a region, the relevant 
lengths would be the widths of the grooves and ridges. 

The local pressures shown in Figure 1 1 are a generalization of those in Figure 10 to allow for 
local pressure jumps at the inlet to each groove and ridge. Unlike Case a, the upper and 
lower envelopes of the local pressure curves will not coalesce as the number of grooves 
becomes large since the magnitudes of the groove and ridge widths do not enter directly into 
the jump relationships given by Equations (3-16) - (3-18) and the magnitudes of the jumps 
will not approach 0 as long as Case b conditions prevail. The local pressure jumps for this 
case will be shown to give rise to a term of order (p*R*N g ) which will be retained in the 
analysis. As with case 1, the other inertia terms will be neglected in the local analysis but 
maintained in the global analysis. These assumptions will thus give rise to asymptotic 
relationships as the number of grooves becomes very large in contrast to the limiting 
relationships corresponding to the assumptions used in Case a. 

c. Jet extends over entire groove. As the number of grooves becomes very large a 
limiting flow should in fact occur as a jet forms that extends across the entire groove. 

The bulk flow theory which is based on the assumptions that the shear forces on the 
turbulent core arise directly from interactions between the core and the bounding surfaces 
break down at this point. A rigorous analysis would involve consideration of the jet, the 
wake, interactions between jet and wake, interactions between jet and wall etc. and 
would require three dimensional solutions to the momentum and continuity equations 
which is beyond the scope of the present effort. An alternative approach would involve 
the use of additional empirical relationships to characterize these interactions so that 
limiting behavior can be obtained. Examples of the use of such relationships are given in 
Reference 23 for labyrinth seals and in Reference 24 for straight through screw seals. 

The narrow groove analysis contained in the rest of this section will refer to Case b unless 
otherwise specified. The theory requires the development of expressions for the local quantities 
in terms ol continuous global quantities that become valid as the number of grooves becomes 
large. Local quantities will be denoted by the subscript g when referring to grooves and r when 
referring to ridges. Variables such as p,u,v appearing without r or g subscripts will henceforth 
denote global quantities. The subscript r denoting ridges has been used for consistency with 
existing literature and should not be confused with the radial position variable r which will not 
be used as a subscript. 

Spiral groove patterns for a given region are characterized by three parameters which are defined 
in terms of quantities shown in Figure 1 1 . These are the groove depth 6 = h„ - h r , the groove to 
pitch ratio a = A0 g /(A0 g +A0 r ) and the groove angle P detined as the angle between the groove 
and the surface velocity direction, measured from the groove. 
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One may now introduce the dimensionless groove and ridge flows in the 0 and s directions 


■g 9 


u h 
g g 


W 1 . 
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5 , h t 


3-21 


and express the turbulent shear functions <5 and Y defined by Equations (3-14) and (3-15) for 
the grooves as 


0 g . - qj(3?i,5si,h a ) 


h„ h 


3-22 


and for the ridges as 


, l|J r - tU(il,^,h r ) 

hf h r hj. hp 


3-23 


One may now write Equations (3-1 1) and (3-12) for the groove regions as 
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and for the ridge regions as 
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The unit tangent and normal vectors for the logarithmic spiral as shown in Figure 1 1 are given 
by 


tp - cos p i + sin (3 j , n p - sin p i - cos p j 


3-28 
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If local inertia effects were neglected, the pressure would be continuous at each groove-ridge 

interface and the derivative of the pressure in the direction tangential to the interface, Vp* t p , 
would also be continuous. In keeping with the assumption that the only local inertia effect to be 
considered is the pressure jump normal to the interface, continuity of the tangential pressure 
derivative will be maintained even when this effect is considered. Continuity of the tangential 
derivative of pressure at the interface may be expressed as 

C0S P Pg.e + sln P Pg, s - C0S P P r ,e + sin P P r ,s ■ 3-29 

One may now substitute Equations (3-24) - (3-27) for the corresponding pressure derivatives in 
the above equation to obtain 


cosp <t> g + sinp - cosp 0 r + sinp tP r 


3-30 


The flow rate relative to a groove-ridge interface in the direction normal to that interface, q n > 
must be continuous. For grooves on the stationary surface (b), that flow rate would be hu*n p 

and for grooves on the moving surface (a), it would be h(u - rco 1 )• n p . One may use the 
dimensionless groove and ridge flow rates defined by Equation (3-21) and perform the indicated 
dot products to express this continuity condition as 

9 n - (q g 0-™h g ljsinp - q gs cosp = (q^ -fwh r IJsinp - q re cosp , 3-31 

where I u is equal to 1 if the grooves are on the moving surface and 0 if the grooves are on the 
stationary surface. The two equalities on the right may be rearranged to obtain the relationship 


sinp q gQ - cosPq gs - sinp q r0 + cospq rs - ftB6l u sinp . 3-32 

One may now equate flow the across a groove-ridge pair in the 0 direction to the global flow in 
that direction and express the result in dimensionless form as 

a V * (1 -°0q r8 - Uh . 3-33 

A similar balance in the s direction results in the relationship 

a Ags + 0 - a )q rs “ vh . 3-34 

The global film thickness, h is the value which results in the same flow area as that produced by 
the local film thickness variation and may be expressed in dimensionless form as 

h - ah g . (1 - a)h r - h r * a6 3.35 
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Equations (3-30) and (3-32) - (3-34) represent four algebraic equations which may be solved for 
the four local flow rates q g0 , q gs , q r0 and q rs . Iterative solutions are obtained numerically with the 
use of a Newton-Raphson procedure to linearize Equation (3-30). 

The pressure changes (downstream - upstream) incurred at the entrance to a groove, Ap 4 „ and at 
the entrance to a ridge, Ap r , may be expressed in dimensionless form with the aid of the jump 
function x, defined by Equations (3-16) - (3-18), as follows: 


- |pR'X(ph fl ,S) , Ap r 


pR-x(^,iv-5) 

h. 


3-36 


The normal flow component, q n , may be expressed in terms of the global velocities and film 
thickness with the use of Equations (3-31) and (3-33) - (3-35) as 


q n - (Gh-rtohlJsinP - vhcosp . 3-37 

By equating the change in pressure over a groove-ridge pair to the changes in the global 
pressure, first in the 6 direction, as shown in Figure 11, and then in the s direction and 
expressing the result in dimensionless form, the following equations are obtained: 


lap 

r ae 


a Pfl,8 + O -°0P,9 + 


AiyAPr 

rA0 


sign(q n ry7) 


3-38 


ap 

as 


a P g .s * (1-a)p r , s 


A FV A Pr 

fA0 1 tan p| 


sign(q n n p j) 


3-39 


where the sign functions appearing in the above equations are either 1 or -1 depending on the 
direction of the normal component of flow with respect to the corresponding coordinate direction. 

Finally, noting that A0 = 2Tc/N g and the definition of n p given by Equation (3-28), one may 
substitute Equations (3-24) - (3-27) for the corresponding local pressure gradients and Equation 
(3-34) for the pressure jumps in Equations (3-38) and (3-39), and write the resulting equations in 
the form 


1_ 6p 

f a0 


pf0(u,v,h,I f ) 


ao + u au N 

^ a t ds r a0 j 


3-40 


ap 

as 


pr4J(u,v,h,l f ) + pfR 




as 


U dV 

f a0^ 
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where 


0 - a$ g * (1 -a)d> r + R * \\~ ^sign(q n sinP) 


N 


f Am 


X(^,h g ,5) * X (?.h rI -5) 


% 


L o 


and 


3-42 


ip . a4J g + (1 -a)ip r * R*fl f ^L 


N 


r 4nr | tan p 


sign(q n cosP) 


X(f,h„,6) . x(?.h,.-6) 
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The Coriolis and centrifugal inertia terms for face seals and the spiral groove local inertia term have 
been combined with the turbulent shear resistance terms in Equations (3-42) and (3-43). This has 
been done for convenience because of their algebraic nature and will simplify the development and 
solution of the first and second order flow equations to be discussed in the next section. 

The actual number of grooves, N g , does not usually appear in a narrow groove analysis. It enters 
here only to characterize the cumulative effects of losses due to sudden changes in cross section 
between ridges and grooves. Formally setting N g =0 would therefore only eliminate the 
inclusion of these local inertia effects and would maintain all of the other spiral groove effects 
under consideration. Caution should be exercised when groove or ridge widths become the same 
order of magnitude as the clearance. This would invalidate one of the basic geometric 
assumptions of the bulk flow theory used throughout this analysis and could have a particularly 
pronounced affect on the validity of the above local inertia terms. 

Equations (3-40) and (3-41) represent the generalization of Equations (3-11) and (3-12) to spiral 
groove seals and reduce to Equations (3-11) and (3-12) when both a and N g are set equal to 0. 
The continuity equation [Equation (3-13)], and the continuity and boundary conditions given by 

Equations (3-16) - (3-20) remain unchanged subject to the interpretation of u, v, p and h as their 
global values when a > 0. 


Development of First- and Second-Order Equations for Small Displacements 

The first-order variables will be the global film thickness, pressure and tangential and transverse 
velocities for an undisplaced shaft. They will depend only on s and their dimensionless forms 

will be denoted by H(S), P(S), U(S) and V(S) respectively. Substitution of the above variables 

for h, p, u, and v respectively, in Equations (3-40), (3-41), and (3-13) yields the following set of 
ordinary differential equations: 

0(U,V f H,l f ) ♦ RV— - 0 , 3_ 44 


dP 

dS 


pPUiCU.V.H,^) 


pc RV — 
dS 


3-45 


d(rHV) Q 
dS 


3-46 
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The boundary and jump conditions for the first-order variables are prescribed by Equations (3-16) - 

(3-20) with h, p, u, and v replaced by the corresponding first-order variables. The boundary 
conditions are given by 


U - u,. 


Pin 


^PR-X(V,H,-~) 


P in - IjSrFtO.OV 2 @S-S it 
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and 


P'Pex @S-S ex . 


3-48 


One may now add a small displacement, h'(s,0,t) to H(s) and express the resulting film 
thickness, pressure and velocities in dimensionless form as 


h-H(S) + h'(S,0,t) , u - U(S) + u'(S,0,t) . 
v - V(S) + v'(S,0,t) , p-P(S).p'(S,0,t) . 


3-49 


When above relationships are substituted for the corresponding variables in Equations (3-40), (3- 
41) and (3-13) and quadratic terms in second-order (primed) quantities are neglected, the 
following second-order, linearized momentum and continuity equations are obtained: 


1 3p' 

- R 

' du' 

. V M ' 

dU - 

+ V 

fpr 30 


K at 

as 

dS 

1 dp' 

- R 

3v' 

♦ V 55 ' 

dV-, 

+ — v 

pr 3S 

* 

k at 

3S 

dS 


0 30 ^ 
f 30 

0 dv' ' 

f 30 


+ %u' + ♦ d>ph' , 

+ Wo*' + . 
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3-51 


vK5L> + ^fh' + rH — 
3S dS 3S 


dS 


30 


H — 
30 


** / 
~ r ^L 

3t 
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The subscripts U, V and H are used to denote partial differentiation with respect to those 
variables. 


The continuity requirements for p' and v' at values of s where h is discontinuous are obtained by 
performing a similar linearization on Equation (3-16) and are given below: 


VAh 

H(H-Ah) 


(H-A h) 
H 


vj , Ap' * -p^R*(XvV / + XHh / ) @S - S 
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The tangential velocity u' will be continuous. The end conditions given below are obtained by 

linearizing Equations (3-19) and (3-20) and substituting the values of U, V and P prescribed by 
Equations (3-47) and (3-48): 


u' - 0 , p' . |pfR(XvV'*XHh') @s - s in 

and 

p' - 0 @ s = s ex . 


3-54 
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Solution to First-Order Flow Equations 

The first-order flow equations are solved by a fairly straightforward numerical procedure. To 
begin with one may integrate Equation (3-46) and express the result in the form 


where the dimensionless transverse inlet velocity V in , is unknown and must be initially assumed. 

One may compute V(S), from the above equation, for the assumed value of V in . The tangential 

velocity, U(S), may now be determined by solving Equation (3-44) numerically subject the first 

boundary condition in Equation (3-47). Once U(S) and V(S) have been established, P(S) may be 
determined by numerical integration of Equation (3-45) subject to the second boundary 

condition in Equation (3-47). The correct value of V in will be that which results in a pressure 
distribution that satisfies Equation (3-48) and is determined numerically by Newton's method. 

The transverse flow rate, Q, may then be calculated as Q = 27ir 0 CV 0 r in H in V m . 

A semi-implicit algorithm has been used for solving Equation (3-44). Examples of semi-implicit 
algorithms may be found in many texts on numerical methods such as Reference 25. Although 
they are slightly more complex than some of the widely used explicit methods such as Runge- 
Kutta, they provide increased stability for stiff' systems such as that which may occur here at 
small values of R . Although Equation (3-44) is a single first-order differential equation, the 
algorithm is presented below for a system of N first-order ordinary differential equations, as will 
occur later in the development of the second-order solution. 

The system of differential equations may be represented in the form 


dY, 

dS 


F i (S ( Y 1 ,Y 2 ,... l Y N ) 


i - 1.2, ,N 
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where Y, denotes the dependent variables. If the values of Y; at the new position (S + AS) are 
denoted by the column vector {Y new }, and the values at the old position by { Y}, the equation 
used for computing { Y ncw } may be written as 

{Y new } - (Y} • AS ( [I] - ^[k] ) ’ { F(S.^,Y„Y 2 , ,Y n ) } , 3 _; 

where [I] is the unit diagonal matrix of order N and [k] is a matrix of derivatives whose 
components are 


~F i (S.^,Y 1 ,Y 2 ,...,Y j ,...,Y N ) , i,j - 1,2,. ,N 


3-58 


Solution to Second-Order Flow Equations 

Although the global film thickness H, remains an arbitrary function of the transverse position S, 
the displacement of the center of the seal will be taken to have the form of a small oscillatory 
translation, e x , in the x direction, for a cylindrical seal or e z , in the z direction, for a face seal and 
a small oscillatory rotation about the y axis, i|j, for either type of seal. If the seal oscillates at 

speed Q, the dimensionless film displacement, -h', will take the form 

-h' = [ip 0 Scos6 + l c e xO cos0 + l f e z0 ] cosQt , 3-59 


where 
>c-1 -*T • 


O 


OTo 

v A 


ip-— cosQt , e x 


"xO 


cosQt , e, 
C 


— cosQt . 
C 
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Each of the three terms on the right hand side of Equation (3-59) may be defined as component 
film displacements, 6 k and the equation may be rewritten as 


-h'-£5 k . 

k-1 

The component film displacements may be expressed in complex form as 


3-61 


b k “ 7^( 5 k + 6 k) - 3-62 

where 

5 k - e k n k (S)e i<Jl0lOI> . 3-63 


NAS A/CR— 2004-2 1 3 1 99/VOL3 


56 



The functions 91 (Z) in the above equation and £s(Z) to be used later denote the real and imaginary 
parts of the complex variable Z respectively. The values of e k , r\ k and J k are given below: 


k - 1, 2, 3 

6 k “ <ft). ®x0' ®z0 

n k - s, i cl i f 

V i, i, o 
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Since the second-order flow equations are linear one may superimpose component solutions 

obtained for h' appearing in Equations (3-50) - (3-54) replaced by the various values of -S k . If 
the component solutions for the dimensionless second-order pressure and velocities are denoted 
by Pk> t? and v k respectively, then the values of p\ u' and v' may be expressed as 

P' - EPk • - E Q k ■ v' - £v k , 3-65 

k»1 k-1 k-1 


where 


Pk - 7 »(Pk * Pk) . u k - 4at(Ci k * u k ) , v k - Ist(v k . v k ) 


3-66 


One may now attempt to separate variables by looking for component solutions in the form 


p* - p‘(S)e i(J ‘ 8±ffi > , u* - u k ± (S)e' <J ‘ 9± ™ > , v* - v k 4 (S)e i(J k 0±Ot) 
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When the above expressions for p* u* and v* are substituted for p', u' and v' in Equations (3-50) 
- (3-53) a set of ordinary differential equations is obtained which may be arranged in the form 




V R 


dU 

dS 


A 


v >f + T^Pk*= e k n k d>^ , 

prr 
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dv 


ds 


k Au* * 

P fH dS k 


e k 

v d( A) . 

( - ) 

dV . a± 
— + iO 

ITI. 

fH 

dS 

l dS J 

•k 
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1^1 
p ds 


v] 

iP 0 -iRV 


UV* 

r dV 

V 


r ) 



l dS 

rH 

dS 

/ 


- e k^HH k - e k R 


fH 


dS 


_dV ■ A 
dS 


♦ iQ* 


r n k 
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where 


3-71 


f 

The transformations given in Equation (3-66) satisfy the periodicity requirements in 8. The 
remaining boundary and continuity conditions are obtained by substituting the variables defined 
in Equation (3-67) for the corresponding variables in Equations (3-53) - (3-55). The form of the 
equations remains unchanged and the results are given below: 


v? 


VAh 


H(H-Ah) 


e k>V 


(H - Ah)-± 

- v kJ 

H 


Ap* - IprR-fXvV^-XHejrik) @S - Sj , 


® » Pk “ 2 P'^(XvV| ( -XfiekHk) @S - S,„ 


and 

P* - 0 @ S - S ex . 


3-72 


3-73 


3-74 


Equations (3-68) - (3-70) are solved numerically using the algorithm defined by Equations (3-56) - 
(3-58). Complimentary and particular solutions are constructed and combined so that the boundary 
conditions given by Equations (3-72) - (3-74) can be satisfied without having to iterate. 

The real component pressures p k , can be expressed in terms of the real and imaginary parts of p k 
and j| by expanding the complex pressures, p k , substituting the result for the complex pressures 
in Equation (3-66) and extracting the real part as indicated. The result may be written as 


P k - (Pk + Pk ) cosJkQ - 3(p; + p k )sinJ k 0 cosOt - ^[s(p; - p k )cosJ k 0 - 91 (p; - p k )sinJ k 0 sinQt . 3_75 


The first term on the right in the above equation will be in phase with the displacement and 
related to the stiffness of the seal. The negative of the second term will be in phase with 
oscillatory velocity and related to the damping coefficient. 


Forces, Moments and Dynamic Coefficients 

The dimensionless applied loads and moments obtained by integrating the pressure over the seal 
area are 


2tt s r 2n s r 2n s r 

w x -l c / f pcos 0 fdSd 0 , w y -l C J | psin 0 fdSd 0 , J prdSd 0 , 

0 s L 0 S L 0 s L 


2n s r 

J psin 0 rSdSd 9 , 

0 s L 


2n s r 

M y - J J pcos 0 fSdSd 0 . 

o s L 
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All first order applied loads and moments will be zero with the exception of the first-order axial 
load acting on a face seal, W z , which is given in dimensionless form by 

Sr 

w z = 2nl f j* PfdS . 3_ 77 

s L 

The various stiffnesses can be obtained by substituting the first term on the right hand side of 
Equation (3-75) for the pressure, p, appearing in the force and moment integrals in Equation (3- 
76) and dividing the result by the associated displacement of the center of the seal. The 
dimensionless translations and rotations that have been included in Equation (3-59) are a rotation 
about the y axis of magnitude ^cosOl: (corresponding to k=l), a translation in the x direction 
of magnitude e x0 cosQ t (corresponding to k=2) for a cylindrical seal and a translation in the axial 
direction of magnitude e z0 cosQ7 (corresponding to k=3) for a face seal. The damping 
coefficients may be obtained in a similar manner by substituting the second term on the right 
hand side of Equation (3-75) for die pressure in Equation (3-76) and dividing the result by the 
velocity -Q t $ 0 sinQ t for k=l, -QS x0 sinQ7 for k=2 and -QS z0 sin£H for k=3. 

For small displacements about the centered position, stiffness and damping values associated 
with translations in the y direction or rotations about the x axis may be obtained from symmetry 
considerations. For example, a displacement in the y direction will produce a force in that 
direction of the same magnitude as that which would be produced in the x direction for an 
equivalent displacement in that direction (K yy =K xx ). Similarly, a displacement in the y direction 
will produce a force at 90 from that direction (which would be the negative x direction) as 
would be produced in the positive y direction for an equivalent displacement in the x direction 

(■^xy --^-yx ) 1 

The symbols c|> and ifr, when appearing as second subscripts on the dynamic coefficients, will 
denote rotations about the x and y axes respectively. When appearing as first subscripts they 
will denote moments about the x and y axes. The dimensionless stiffness and damping 
coefficients relating moments about the x and y axes to a rotation about the y axis (k=l) for both 
face and cylindrical seals are 


k * v '2tf S (Pf'*Pi> SdS ■ 9t(P,’.pi)SdS , 

1 S L 1 S L 

S R S R 

/ 8t(Pi P, )SdS , B --H- f 3(p,-p-)SdS 
2Dr. i one J 
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i s, 


The conesponding coefficients for a rotation about the x axis obtained from symmetry 
considerations are 



- K44 ■ - B w , - -B 
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The dimensionless dynamic coefficients relating the x and y force components to a rotation 
about the y axis (k=l) for a cylindrical seal are 


ttI , R _ _ - nl _ _ 

V 2 f f ^(Pi - Pi> dS . K f S(p; + Pi)dS , 

1 S L 1 s L 

3-80 

nl _ _ - nl S . R _ _ 

/ 3( P ;- Pl -)ds , I a( Pl *- Pl -)ds 

2^ e i s L 20e 1 Sl 


and the corresponding coefficients for a rotation about the x axis are 


^x4> " » ^y$ ' ^XIJJ > ^X(J> “ ^yl)J * ^y$ “ ® x ij» ' 3-81 

The dimensionless stiffness and damping coefficients relating moments about the x and y axes to 
a translation in the x direction (k=2) for a cylindrical seal are 

S R S R 

/ 3(p 2 .p 2 )SdS , K f 8t(pj*p 2 )SdS , 

2 S L 2 s L 
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S R S R 

B*-— f St(p 2 -p 2 )SdS , B .-g- f S(p 2 -p 2 )SdS 

2De 2 i L 20e 2 { 

and the corresponding coefficients for a translation in the y direction are 


^Sj>y “ ' ^iyy = K*x ' ®*y “ ®ijjx ’ “ ^4>x ' 3-83 

The dimensionless dynamic coefficients relating the x and y force components to a translation in 
the x direction (k=2) for a cylindrical seal are 


nl _ _ „ nl R 

^oc— /#(Pi*Pi>dS . ^-—2- /3(K*K)dS . 

2 S L 62 S L 
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TT| | Sp 

>*x ■ -4- / S(p 2 -p 2 )ds , B^-4 5 - / «(p 2 -p 2 )dS 
2Qe 2 Sl 20e 2 Sl 


and the corresponding coefficients for a translation in the y direction are 

^xy ' “Kyx ’ ^yy ^xx * ^xy ~^yx ’ ^yy ' ^xx ' 
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Finally, the stiffness and damping coefficients relating the axial force to an axial translation 
(k=3) for a face seal are 


nl R nl R 

Kzz f 9t(P3*P 3 )dS , - - r -!~ [ S(p 3 -p 3 )dS . 

e 3 s L Oe 3 { 

As a result of the assumption of small displacements about the centered position used in this 
analysis, the remaining dynamic coefficients for a face seal will all be zero: 
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K* ' ' Kk ' - B* z - B zv - - 0 . 3-87 

Relationships of the following type may be used to calculate the physical stiffness and damping 
coefficients from the dimensionless ones 


If Poh) |/- Po r O 1/ \y 

Nx ■ Nx . Ah > •Sc* 


Po r o 
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and 
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- r 3 
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B, 
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B„ 


P° r o § 
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X<t> p\/ “x* 

0 v -' v o L ' v o 

The stiffness coefficients, which relate to forces in phase with the displacements and the 
damping coefficients, jvhich relate to forces 90° out of phase with the displacements, are both, in 
general, functions of Q. For the case of incompressible flow under consideration, the damping 
coefficients and die cross coupled stiffness coefficients (K xy , K^, K yx , K^) are nearly 
independent of Q. Under these circumstances, it has been found that the frequency dependence 
may be described adequately with the introduction of an apparent mass matrix [A] such that 
[K] = [K°] - Q 2 [A], where [K°] is the stiffness matrix evaluated at Q=0. The components of [A] 
may be evaluated by computing [K] at a characteristic speed such as Q=Vo/r 0 corresponding to 
Q=1 or synchronous speed (Q=g)) and using the formula 



3-90 


The cross coupled components of [A] which have been included for completeness are not 
significant and the direct components are indeed found to be nearly independent of Q . 

Although the linearizations that have been used in obtaining the second-order flow equations are 
only valid for small values of e k , once they have been performed the resulting equations for the 
dynamic coefficients may be shown, by rescaling variables, to be independent of the value e k 
actually used. Any convenient value such as e k =l may be used in computations. 


NAS A/CR— 2004-2 1 3 1 99/VOL3 


61 



Determination of Power Loss and Flow 

When the rotating surface (a) is ungrooved, the power loss will be proportional to the integral of 
r v i over the surface of the seal. When the rotating surface is grooved, an additional term will 
be needed to correct for the pressures acting edges of the grooves. Although the solution to the 
equations based on bulk flow theory used here does provide for the prediction the shear stress on 
each surface, the validity of power loss predictions is somewhat doubtful. This is due to the 
difference in characterizing the bulk flow behavior for velocity distributions associated with 
Couette and Poiseuille flow. 

One may write x a as x a = t c + x p where t c denotes the Couette part of the shear stress defined 
as t c = (' r a + x b )/2 and t p denotes the Poiseuille portion of the shear stress defined as x p = (x a - 
t b )/2. It can be seen from Equations (3-1) and (3-2) that only f p appears in the primary and 
secondary flow equations used to predict load, flow and dynamic coefficients. Comparisons 
between theory and experiments based on measurements of flow and dynamic coefficients alone 
may be used to validate or determine values of i p , but will not necessarily result in a bulk flow 
model that accurately predicts power loss. This weakness in the bulk flow model was 
recognized by Hirs (Reference 14) who suggested the use of different friction factor 
representations based on measurements obtained from experiments based on Couette and 
Poiseuille flow. The method for predicting power loss used here will be a highly approximate 
one based on generalizations of an approach used by Hirs for cases where both surfaces were 
smooth and the flow was predominantly in the tangential direction. 

If the scaler shear stress, t, is used to represent the tangential shear stress, (t=t-~i ), then the 
dimensionless forms of the tangential shear stresses x a and x b , can be obtained from Equations 
(3-4), (3-5), (3-8), and (3-10) as 


^ ■ -pRA( R b )“ • 3 . 

h h 

One may attempt to compensate for the influence of the Couette portion of the shear stress by 
replacing t a with the effective shear stress, t*, defined as 


T - 




3-92 


in the integral used for determining power loss. The correct value of the "Couette reduction 
factor", A, may be shown to be A=3 for laminar flow. A value of A=1.2 will be used for 
turbulent flow. This value was suggested by Hirs based on one dimensional Couette flow 
measurements and provides power predictions that are in good agreement with those obtained 
with Ng and Pan theory (Reference 26) over its range of validity. 


The dimensionless power loss, 9 = P/(Cp 0 roO)), obtained from the solution to the first-order flow 
equations is given by 


P - -2" / 

s. 


at* ♦ (1 -a)f* r ♦ l w 6 


aP 


N &R* 
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The quantities T g and T r appearing in the above equation denote the values of the effective shear 
stress, t , based on the first-order velocities and film thicknesses for the grooves and ridges 
respectively. The effective shear stress is defined by Equations (3-91) and (3-92) and the local 
velocities are obtained from the solution to Equations (3-30) and (3-32) - (3-34) and the 
definitions provided by Equation (3-21). The sum of the first two terms appearing inside the 
brackets of the integrand in Equation (3-93) represent the global value of t*. The third term has 
been added to include the forces acting on the edges of the grooves when the grooves are on the 

rotor. The quantity P g>0 is obtained by substituting first-order variables for the corresponding 
quantities in Equation (3-24) and is given below: 


fl.9 


-pfO> - prR 
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3.2 Overview of Computer Code SPIRALI 

SPIRALI has been written to implement the analysis that is contained in the preceding section. 
Full input and output descriptions along with a number of sample runs are given in Reference 
28. A brief discussion of a few items not covered in Section 2.1 is given below. 

Although SPIRALI is designed to run over a full range of conditions of laminar or turbulent 
flow, certain precautions need to be taken when the transverse flow rate is small. It can be seen, 
for example, that Equation (3-44)J)ecomes singular as R*V approaches 0. Numerical problems 
can be avoided here by setting R*V=0 and numerically solving Equation (3-44) as an algebraic 
equation rather than as a differential equation. SPIRALI provides an option for the user to 
exclude the transverse inertia effects arising from the terms vdu/dS and I f uv/r in Equation (3-11) 
and vdv/dS in Equation (3-12) and the associated discontinuities in pressure that occur at the 
inlet and at jump discontinuities in film thickness. This option, which is controlled by the input 
parameter NOI (see input description) is particularly useful for situations where the rotation 
speed is high but the seal pressure difference is low. 

If transverse inertia is included, the location of the upstream or inlet boundary must be specified. 
When the seal pressure difference is sufficiently high, the upstream boundary will be at the high 
pressure side of the seal. At lower pressure differences the upstream boundary may depend on 
the rotor speed and the spiral groove parameters. The input parameter IFLOW is used to select 
the upstream boundary. If one is not sure of the flow direction, he may first run the program 
with transverse inertia excluded. One may then include transverse inertia by adjusting IFLOW 
in accordance with the sign of the transverse flow rate given in the preceding output and 
rerunning the program with transverse inertia included. 

The analytical procedure contained in Section 2 has been oriented toward determining pressure 
distribution, load, flow, power loss, stiffness and damping for a given centered film thickness 
distribution. Face seals will, in general, have a non-zero first-order axial load component and it 
is often desirable to determine the equilibrium film thickness from a prescribed load used to 
balance the seal. SPIRALI provides an iterative homing option for this purpose that is 
controlled by the input parameter IHOME. The homing option should be used with caution 
since convergence is sensitive to the starting film thickness which is determined from the input 
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parameter C. It is recommended that the program be run for various values of C without the homing 
option so that a good starting value of C (one that produces an output value of the load to balance the 
face seal that is reasonably close to the prescribed load) can be obtained for the homing process. 

Although the analysis provides for an arbitrary transverse film thickness variation, SPIRALI is 
limited to a quadratic film variation with coefficients determined from the parameters H TAP , H BRL 
and C as shown in Figure 12. These parameters may be positive, as shown or negative and 
should be able to characterize a fairly general film variation without the complexity of required 
inputs at every grid point. 

The logic used in SPIRALI for solving the first and second-order flow equations, computing 
load, flow, power loss, stiffness, damping, apparent mass coefficients, etc. for a given case is 
shown in Figure 13. 

3.3 Verification 

An available MTI computer code, PN471, has been used to obtain a partial check on the results 
of SPIRALI. This code implements a two dimensional finite difference solution to the spiral 
groove equations that is quite different from the perturbation methods used in SPIRALI. It does 
not treat inertia effects and uses the Ng and Pan (Ref. 25) linearized turbulence theory which 
limits its validity to cases where turbulence is dominated by circumferential component of the 
flow. The assumptions used in the MTI code will be identical to those used in SPIRALI when 
the flow is laminar and inertia effects are neglected. The bulk flow model used in SPIRALI 
should also produce similar (although not identical) results to those obtained with the MTI code 
under turbulent conditions when inertia effects are neglected, there is no significant seal pressure 
difference, the circumferential Reynolds numbers are between 10 4 and 10 5 and the transverse 
Reynolds numbers are small compared with the circumferential ones. 

Four test cases obtained with SPIRALI are presented for comparison with the MTI code. Cases 
1 and 3 correspond to laminar flow in a spiral grooved cylindrical and face seal respectively. It 
can be seen that the results of both cases agree extremely well with the results of the MTI code. 
The small discrepancies can easily be explained by truncation error associated with the finite 
difference solution used in the MTI code. Cases 2 and 4 provide comparisons similar to cases 1 
and 3 under conditions of turbulent flow. It can be seen that the Reynolds numbers are in the general 
range of validity of the MTI code and all results agree to within 10% for both cases which is 
acceptable considering the differences in the turbulent flow models. 

While the above cases provide a partial check on turbulent, spiral groove seals, they do not 
provide any validation on the treatment of inertia effects or seal behavior under high imposed 
pressure gradients. In order to obtain verification of the treatment of these effects, comparisons 
are presented between the results published by Childs in Reference 17 and those obtained with 
SPIRALI for a plain ungrooved cylindrical seal. The output from SPIRALI is given as cases 5-8 
and comparisons of flow rates and rotordynamic coefficients with those published by Childs are 
given in Table 1. Although agreement is fairly good, discrepancies in results of up to 10%, such 
as that which occurs in the direct stiffness for case 6 where u in =r 0 oo/2 and L/D=l, are somewhat 
surprising since the assumptions for that case are identical with those used by Childs. In order to 


♦SPIRALI test cases and output are presented as cases 1-8 beginning on page 78. 
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Figure 13. Mow Diagram for Overall Logic Used in Computations 
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obtain a further independent check on the results, a reduced form of the differential equations 
published by Childs for predicting the direct and cross stiffness coefficients when u in =r 0 o)/2, were 
programmed and solved by the Runge-Kutta method. The program produced values of K xx , K yy and 
Q that agree to four places with the case 6 results produced by SPIRALI. It is thus likely that the 
small discrepancies between results obtained from SPIRALI and those published in Reference 17 are 
due to truncation errors in the latter solution. 

Parallel groove geometries can be analyzed by using the separated region option of SPIRALI. 
Analysis was conducted of the parallel groove and stator 1 of Reference 27 by treating each groove 
and ridge as a separate region. The clearance and groove geometry are shown in Figure 14. 

A total of 19 regions were used consisting of nine grooves, eight lands between grooves, an inlet 
land and an exit land. (The parameter statements in SPIRALI governing the maximum number of 
regions were enlarged appropriately). The inlet and exit land lengths were enlarged slightly to make 
the sum of the lengths of the individual regions equal to the total seal length (50.8 mm). The seal 
radius is 50.8 mm. The viscosity and density are 1.54 x 104 Pa-s and 1570 kg/m 3 , respectively. The 
program default values for flow parameters of m = -0.25, n = 0.0791 (Hirs coefficients for smooth 
surfaces as given by Blasius theory), and loss coefficient of ( = 0 were used. Although it was stated 
in Reference 27 that no intentional swirl was provided, a swirl speed of 25% of the rotor speed was 
used in all computations to characterize the swirl that might naturally be present. 


Comparison of Flow Rates and Force Coefficients 

Experimental flow rates are presented in terms of a discharge coefficient, C D , defined by 
Equation (3) of Reference 27. The value of C D 1/2 for stator 1 is shown in Figure 6 of Reference 
27 to have a constant value of approximately 0.5 for all measured pressures. A comparison 
between that value and those predicted by SPIRALI is shown in Figure 15. The agreement is 
good and the results are predicted to be relatively insensitive to the rotating speed over the range 
of conditions covered, as indicated in Reference 27. 

The dynamic coefficients presented in Table 1 of Reference 27 are computed from measured force 
coefficients based on the assumption that the coefficients are independent of speed. For purposes of 
comparison, the measured force coefficients, f 9 and f r , were retrieved with the use of Equation 6 in 
Reference 27. The seal pressure drops, AP, were calculated from the axial Reynolds number, R a , 
using Equation (3) of Reference 27 with C D 1/2 = 0.5. Values of f e versus AP obtained in this manner 
and those predicted by SPIRALI appear to be in good agreement, as shown in Figure 16. The 
predicted variation in f r with rotating speed (Figure 17) is somewhat weaker than that obtained by 
experiment. Discrepancies in f r can, in part, be due to sensitivity to the input values to SPIRALI 
provided for the friction factor and loss coefficients. Also, the effective mass, M d , in Table 1 of 
Reference 27 shows a significant jump at a Reynolds number of 200,000 that may be suspect. 

Effective stiffness values are compared in Figure 18. The experimental stiffnesses are based on 
extrapolations to zero speed. Stiffness is assumed to be caused principally by the pressure gradient 
and inlet Bernoulli effect which causes a circumferential pressure gradient and positive stiffness if 
the rotor goes eccentric. As Figure 18 indicates, the predicted values approach the experimental 
values as the speed is reduced. The predicted decrease in K with increasing rotating speed is to be 
expected due to the Bernoulli effect associated with fluid rotation. A corresponding comparison of 
C cf with C - k/ co is given in Figure 19. In both cases, the agreement improves as the rotational speed 
is decreased and the physical interpretation of the data becomes valid. 
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Figure 14. Parallel-Groove Pressure Breakdown Seal 



AP (bar) 


I'igure 15. Parallel-Groove Seal Flow Coefficient 
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Figure 1 7. Normal Force Coefficients 
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Figure 18. Parallel-Groove Seal Effective Stiffness 
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Figure 19. Comparison Between C - k/co and C ef at Various Rotating Speeds 
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Effect of Local Pressure Discontinuities on Predicted Flow 

As described above, the effects of pressure jumps at the grooves are accounted for by SPIRALI 
in the prediction of flow and force coefficients. An estimate of these terms can be obtained by 
comparing the flow rate for circular grooves with that which would be predicted to occur if the 
effects of pressure jumps at the groove interfaces were neglected. This comparison is made at 
1000 rpm in Figure 20. The large discrepancies between the two curves indicate a significant 
effect of local inertia that must be accounted for. 

Check cases were also run against the helical groove results of Reference 27. Figure 21 shows 
comparisons of the flow coefficient for varying groove angles at a loss coefficient equal to 1, and the 
theory predicts results accurately. Cross-coupled force comparisons are shown on Figure 22. Cross- 
coupled effects are important because they are a measure of the stability of the seal, and in some 
cases, cross-coupled forces cause subsynchronous whirling of the rotor. The comparative results are 
very good. It is noted that cross-coupled forces were compared rather than the stiffnesses indicated 
in Reference 27. The experimental procedure measures forces and computes stiffness and damping 
based on the assumption that the cross-coupled coefficients vary directly with speed. The 
assumption of zero cross coupling at zero speed is suspect for helical grooves because the grooves 
apply a pressure-driven tangential flow at zero speed opposite the direction of rotation. As shown in 
Figure 23, the comparisons of direct stiffness, K ef , were not as good as the parameters presented 
above. This may be due to the experimental assumptions or code assumptions. 

The net results of these studies is that SPIRALI is an effective tool for analyzing many parallel 
and helical groove seals. The principal limitation is the bulk flow model that does not treat flow 
variations across gaps. If large gaps are used, then more sophisticated CFD analysis is required. 



Figure 20. Effect of Local Pressure Discontinuities on Predicted Axial Flow Rates 
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Figure 21. Flow Coefficient : SP1RALI Compared to Reference 27 
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( CASE 1 ) Cylindrical seal with grooves, laminar, no press, grad. 
CYLINDRICAL SEAL, INERTIA NEGLECTED 

LENGTH, DIAMETER, CLEARANCE = 5.0000E-01, 2.0000E+00, 1.0000E-03 (IN) 

ROTOR, SWIRL AND DIST . SPEEDS = 5.0000E+04, 2.5000E+04, O.OOOOE+OO (RPM) 

PRESSURE AT START, END AXIAL BOUNDARIES - 0.0000E+00, O.OOOOE+OO (PSI) 

VISCOSITY = 3.0000E-08 (PSI-SEC), DENSITY = O.OOOOE+OO ( LB - SEC 2 / IN 4 ) 

FLOW = 1.1506E + 00 (IN 3 /SEC) 

TORQUE = 4.3909E-01 (IN-LB), FILM POWER LOSS = 3.4835E-01 (HP) 

AXIAL REYNOLDS NUMBER = O.OOOOE+OO 

CIRC. REYNOLDS NUMBERS FOR ROTOR AT SEAL ENDS - O.OOOOE+OO, O.OOOOE+OO 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP. 

x ( IN) 

y (IN) 

p (RAD) 

p (RAD) 

. FORCE UNIT 

Kx 

2 . 7021E+04 

1 . 3511E + 04 

-6 . 3081E+02 

-1 . 2959E + 03 

LB 

Ky 

-1 . 3 511E + 04 

2 . 7021E + 04 

1 . 2959E + 03 

-6 . 3081E+02 

LB 

Kp 

2 . 3387E+02 

- 6 . 5242E + 01 

1 . 6699E+02 

8 . 6419E + 01 

IN-LB 

Ki|r 

6 . 5242E+01 

2 . 3387E + 02 

-8 . 6419E+01 

1 . 6699E + 02 

IN- LB 

Bx 

5 . 12 97E+00 

7 . 8477E-11 

-2 . 0242E-02 

1 . 3973E-01 

LB-SEC 

By 

-7 . 8477E-11 

5 . 1297E + 0Q 

-1 . 3973E-01 

-2 . 0242E- 02 

LB-SEC 

Bp 

-2 . 0242E- 02 

-1 . 3973E-01 

3 . 0437E-02 

3 . 8319E- 13 

IN-LB-SEC 

Bp- 

1 . 3973E-01 

-2 . 0242E-02 

-3 . 8319E-13 

3 . 0437E- 02 

IN-LB-SEC 

Ax 

4 . 7708E-19 

4 . 0552E- 18 

7 . 4544E-20 

-1 . 7891E-19 

LB-SEC 2 

Ay 

-4 . 0552E-18 

4 . 7708E-19 

1 . 7891E-19 

7 . 4544E-20 

LB-SEC 2 

Ap 

1 . 1927E-19 

5 . 4045E-20 

-1 . 8636E-20 

3 . 3545E-20 

IN-LB-SEC 2 

Ap 

-5 . 4045E-20 

1 . 1927E-19 

-3 . 3545E-20 

-1 . 8636E-20 

IN-LB-SEC 2 


RESULTS OF MTI SPIRAL GROOVE CODE FOR COMP AR I S ON WITH CASE 1 
FLOW = 1 . 1506E + 00 (IN 3 /SEC) 


TORQUE = 4 . 3909E-01 (IN-LB), FILM POWER LOSS = 3.4835E-01 (HP) 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP . 

x (IN) 

y (IN) 

p (RAD) 

p (RAD) 

FORCE UNIT 

Kx 

2 . 6 996E + 04 

1 . 3497E+04 

-6 . 2 96 6E + 02 

-1 . 2947E + 03 

LB 

Ky 

-1 . 3497E+04 

2 . 6 996E + 04 

1 . 2947E + 03 

-6 . 2966E + 02 

LB 

Kp 

2 . 3429E+02 

-6 . 4947E + 01 

1 . 6666E + 02 

8 . 6250E+01 

IN- LB 

Kp 

6 . 4947E+01 

2 . 3429E+02 

-8 . 62 50E + 01 

1 . 6666E + 02 

IN-LB 

Bx 

5 . 1280E+00 

2 . 4507E-07 

-2 . 0178E- 02 

1 . 3967E-01 

LB-SEC 

By 

9 . 3415E-09 

5 . 1280E + 00 

- 1 . 3967E-01 

-2 . 0178E- 02 

LB-SEC 

Bp 

-2 . 0178E-02 

-1 . 3967E-01 

3 . 0411E- 02 

3 . 0462E- 09 

IN- LB-SEC 

Bp 

1 . 3967E-01 

-2 . 0178E-02 

4 . 9311E-09 

3 . 04 HE- 02 

IN-LB-SEC 
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( CASE 2 ) Cylindrical seal with grooves, turbulent, no press, grad. 
CYLINDRICAL SEAL, INERTIA NEGLECTED 

LENGTH, DIAMETER, CLEARANCE = 5.0000E-01, 2.0000E+00, 1.0000E-03 (IN) 

ROTOR, SWIRL AND DIST. SPEEDS = 5.0000E+04, 2.5000E+04, 0.0000E+00 (RPM) 

PRESSURE AT START, END AXIAL BOUNDARIES = O.OOOOE+OO, O.OOOOE+OO (PSI) 
VISCOSITY = 3.0000E-08 (PSI-SEC) , DENSITY = 1.0000E-04 (LB-SEC 2 /lN 4 ) 

FLOW = 1 . 7348E+00 (IN 3 /SEC) 

TORQUE = 6 . 13 00E+00 (IN-LB), FILM POWER LOSS = 4 . 8631E+00 (HP) 

AXIAL REYNOLDS NUMBER - 1.8407E+03 

CIRC. REYNOLDS NUMBERS FOR ROTOR AT SEAL ENDS = 3.8561)3+04, 1.7453E+04 


DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP. 

x (IN) 

y (IN) 

0 (RAD) 

d (RAD) 

FORCE UNIT 

Kx 

2 . 0689E+05 

1. 0073E+05 

-1 . 9681E+03 

-9 . 2019E+03 

LB 

Ky 

-1 . 0 073E+05 

2 . 0689E+05 

9 . 2019E+03 

-1 . 9681E+03 

LB 

K4> 

8 . 8445E+02 

-2 . 4291E+02 

1 . 1123E+03 

4 . 9214E+02 

IN-LB 

KlR 

2 . 42 91E+02 

8 . 8445E+02 

-4 . 9214E+02 

1 . 1123E+03 

IN-LB 

Bx 

3 . 7161E+01 

3 . 3996E- 07 

-6 . 9734E- 02 

5 . 2914E-01 

LB-SEC 

By 

-3 . 3996E-07 

3 . 7161E+01 

-5 . 2914E-01 

-6 . 9734E-02 

LB-SEC 

Bd 

-6 . 9734E-02 

-5 . 2914E-01 

1 . 7825E- 01 

2 . 1459E- 12 

IN-LB-SEC 

Bijr 

5 . 2914E-01 

- 6 . 9734E- 02 

-2 . 1459E-12 

1. 7825E-01 

IN-LB-SEC 

Ax 

3 . 8167E-18 

7 . 6333E-18 

2 . 3854E- 18 

-1. 6698E-18 

LB -SEC 2 

Ay 

-7 . 6333E- 18 

3 . 8167E-18 

1 . 6698E-18 

2 . 3854E-18 

LB-SEC 2 

Ac}) 

1 . 2076E-18 

-3 . 5781E-19 

1 . 4909E-19 

-4 . 4727E-20 

IN-LB-SEC 2 

Af 

3 . 5781E-19 

1 . 2076E-18 

4 . 4727E-20 

1 . 4909E-19 

IN- LB-SEC 2 


RESULTS OF MTI SPIRAL GROOVE CODE FOR COMPARISON WITH CASE 2 


FLOW - 

1 . 7 076E+0 0 

(IN 3 /SEC) 




TORQUE 

- 5 . 8953E+00 

(IN-LB) , 

FILM POWER 

LOSS = 4.677 0E+ 0 0 (HP) 


DYNAMIC COEFFICIENTS ( 

FORCE UNIT / 

DISP. UNIT ) 


DISP. 

x (IN) 

y (IN) 

d (RAD) 

d (RAD) 

FORCE UNIT 

Kx 

2 . 1565E + 05 

1 . 0203E+05 

-2 . 0751E+03 

-9 . 8637E+03 

LB 

Ky 

-1 . 0203E + 05 

2 . 1565E+05 

9 . 8637E + 03 

-2 . 0751E+03 

LB 

Kd 

9 . 7231E+02 

- 2 . 4 006E + 02 

1 . 1627E+03 

5 . 0513E+02 

IN-LB 

Kd 

2 . 4006E + 02 

9 . 7231E+02 

-5 . 0513E + 02 

1 . 1627E+03 

IN-LB 

Bx 

3 . 7602E+01 

2 . 8881E-06 

-6 . 8925E- 02 

5 . 6836E- 01 

LB-SEC 

By 

3 . 6843E-07 

3 . 7602E+01 

-5 . 6836E- 01 

-6 . 8930E-02 

LB-SEC 

Bd 

-6 . 8929E-02 

-5 . 6836E-01 

1 . 8265E- 01 

6 . 4583E- 08 

IN-LB-SEC 

Bi|i 

5 . 6836E-01 

-6 . 8929E-02 

-3 . 9709E-09 

1 . 8265E-01 

IN-LB-SEC 
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( CASE 3 ) Face seal with grooves, laminar, no press, grad. 

FACE SEAL, INERTIA NEGLECTED 

ID, OD, NOMINAL FILM THICKNESS = 1.0000E + 00, 2. OOOOE+OO, 1.0000E-03 (IN) 

ROTOR, SWIRL AND DIST. SPEEDS = 5.0000E+04, 2.5000E+04, O.OOOOE+OO (RPM) 

INSIDE, OUTSIDE PRESSURE = O.OOOOE+OO, 0.0000E+00 (PSI) 

VISCOSITY = 3.0000E-08 (PSI-SEC) , DENSITY - O.OOOOE+OO (LB-SEC 2 /IN 4 ) 

AXIAL LOAD TO BALANCE FACE SEAL = 1.0391E + 01 (LB) 

FLOW = 5.7837E-01 (IN 3 /SEC) 

TORQUE = 2 . 1724E- 01 (IN-LB), FILM POWER LOSS = 1.7234E-01 (HP) 

RADIAL REYNOLDS NUMBER AT ID, OD = O.OOOOE+OO, O.OOOOE+OO 

CIRC. REYNOLDS NUMBERS FOR ROTOR AT ID, OD = O.OOOOE+OO, 0.0000E + 00 


DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) k 


DISP. 

z (IN) 

d (RAD) 

1| / (RAD) 

FORCE UNIT 

Kz 

2 . 3640E + 04 

0 . 0000E+00 

0 . 0000E+00 

LB 

K(f> 

0 . 0000E+00 

6 . 2646E+03 

7 . 0922E+03 

IN- LB 

Kl]/ 

0 . OOOOE+OO 

- 7 . 0922E+03 

6 . 2646E+03 

IN-LB 

Bz 

8 . 9011E+00 

0 . 0000E + 00 

0 . OOOOE+OO 

LB-SEC 

B4> 

0 . OOOOE + OO 

2 . 5921E+00 

2 . 6323E-05 

IN-LB-SEC 

BtJ/ 

0 . OOOOE+OO 

-2 . 6323E-05 

2 . 5921E+00 

IN-LB-SEC 

Az 

0 . OOOOE+OO 

0 . OOOOE+OO 

0 . OOOOE+OO 

LB-SEC 2 

A4> 

0 . 0000E+00 

-1 .1927E-19 

1 . 1927E-19 

IN-LB-SEC 2 

a d 

0 . OOOOE+OO 

-1 . 1927E-19 

-1. 1927E-19 

IN-LB-SEC 2 


RESULTS OF MTI SPIRAL GROOVE CODE FOR COMPARISON WITH CASE 3 


AXIAL LOAD TO BALANCE FACE SEAL - 1.0391E + 01 (LB) 


FLOW - 

5 . 7837E-01 

(IN 3 /SEC) 



TORQUE 

= 2 . 1724E-01 

(IN-LB) , 

FILM POWER 

LOSS = 1.72: 


DYNAMIC COEFFICIENTS ( 

FORCE UNIT / 

DISP. UNIT 

DISP . 

z (IN) 

d (RAD) 

d (RAD) 

FORCE UNIT 

Kz 

2 . 3639E + 04 

7 . 1718E- 03 

7 . 1755E- 03 

LB 

Kd 

5 . 8063E-06 

6 . 2581E+03 

7 . 0834E + 03 

IN-LB 

Kd 

1 . 1824E- 05 

-7 . 0834E+03 

6 . 2581E+03 

IN- LB 

Bz 

8 . 9022E+00 

-2 . 7999E-07 

-3 . 9768E-07 

LB-SEC 

Bd 

-6 . 6333E-08 

2 . 5909E+00 

2 . 9199E- 06 

IN-LB-SEC 

Bd 

-1 . 9 18 IE- 0 8 

-2 . 9131E- 06 

2 . 5909E+00 

IN-LB-SEC 
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( CASE 4 ) Face seal with grooves, turbulent, no press, grad. 

FACE SEAL, INERTIA NEGLECTED 

ID, OD, NOMINAL FILM THICKNESS = 1.0000E+00, 2. OOOOE+OO, 1.0000E-03 (IN) 

ROTOR, SWIRL AND DIST. SPEEDS = 5.0000E+04, 2.5000E+04, 0. OOOOE+OO (RPM) 

INSIDE, OUTSIDE PRESSURE = 0. OOOOE+OO, 0.0000E+00 (PSI) 

VISCOSITY - 3 . 0000E- 08 (PSI-SEC), DENSITY = 1.0000E-04 (LB-SEC 2 /IN 4 ) 

AXIAL LOAD TO BALANCE FACE SEAL - 5.8223E+01 (LB) 

FLOW = 7 . 2459E- 01 (IN 3 /SEC) 

TORQUE = 2.4948E+00 (IN-LB), FILM POWER LOSS = 1/9792E+00 (HP) 

RADIAL REYNOLDS NUMBER AT ID, OD - 1.5376E+03, 7.6881E+02 

CIRC. REYNOLDS NUMBERS FOR ROTOR AT ID, OD = 1.8985E+04, 1.7453E+04 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP. 

z (IN) 

(}) (RAD) 

ijJ (RAD) 

FORCE UNIT 

Kz 

1 . 2740E+05 

0 . 0000E+00 

0 . 0000E+00 

LB 

K4> 

0 . 0000E+00 

3 . 4989E + 04 

3 . 8692E+04 

IN- LB 

Ki|r 

0 . 0000E+00 

-3 . 8692E+04 

3 . 4989E+04 

IN-LB 

Bz 

4 . 8105E+01 

0 . 0000E+00 

0 . 0000E+00 

LB-SEC 

B0 

0 . 0000E+00 

1 . 4225E+01 

6 . 4594E- 05 

IN-LB-SEC 

B\Jr 

0 . 0000E+00 

-6 . 4594E- 05 

1.4225E+01 

IN-LB-SEC 

Az 

0 . 0000E+00 

0 . OOOOE+OO 

0 . OOOOE+OO 

LB -SEC 2 

A0 

0 . OOOOE+OO 

1 . 9083E-18 

-9 . 5417E-19 

IN-LB-SEC 2 

A0 

0 . OOOOE+OO 

9 . 5417E- 19 

1. 9083E-18 

IN-LB-SEC 2 


RESULTS OF MTI SPIRAL GROOVE CODE FOR COMPARISON WITH CASE 4 
AXIAL LOAD TO BALANCE FACE SEAL = 5.8223E+01 (LB) 

FLOW = 7 . 2 032E- 01 (IN 3 /SEC) 

TORQUE = 2 . 4013E+00 (IN-LB), FILM POWER LOSS = 1.9050E+00 (HP) 

DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP . 

z (IN) 

0 (RAD) 

0 (RAD) 

FORCE UNIT 

Kz 

1 . 3685E + 05 

5 . 1323E- 02 

5 . 1174E-02 

LB 

K0 

8 . 1427E-06 

3 . 7351E + 04 

4 . 0251E+04 

IN- LB 

KiJj 

3 . 8 13 IE- 0 5 

-4 . 0251E+04 

3 . 7351E+04 

IN-LB 

Bz 

4 . 9999E+01 

-3 . 3078E- 06 

-3 . 4927E-06 

LB-SEC 

B0 

9 . 6839E-07 

1 . 4768E+01 

7 . 7653E-06 

IN-LB-SEC 

B0 

2 . 8875E-07 

-6 . 9217E- 06 

1 .4768E + 01 

IN-LB-SEC 
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( CASE 5 ) Childs finite length solution, RPM0=RPM/2, L/D=.2 
CYLINDRICAL SEAL 

LENGTH, DIAMETER, CLEARANCE = 3.0480E-02, 1.5240E-01, 1.9050E-04 (m) 

ROTOR, SWIRL AND DIST. SPEEDS = 3.6000E+03, 1.8000E+03, O.OOOOE+OO (RPM) 

PRESSURE AT START, END AXIAL BOUNDARIES = 3.4400E+06, 0.0000E+00 (Pa) 

VISCOSITY = 1 . 2950E-03 (Pa-SEC) , DENSITY = 1.0000E+03 (Kg/m3) 

FLOW - 4 . 0061E-03 (nd/SEC) 

TORQUE = 2 . 2528E+00 (N-m) , FILM POWER LOSS = 8.4929E+02 (WATT) 


DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP . 

x (m) 

y (m) 

<J> (RAD) 

ijr (RAD) 

FORCE UNIT 

Rx 

1 . 8896E+07 

4 . 1269E+06 

-3 . 3389E+04 

1 .4185E+06 

N 

Ky 

-4 . 1269E+06 

1 . 8896E+07 

-1 . 4185E+06 

-3 . 3389E+04 

N 

Kcj) 

-1.2847E+04 

9 . 8025E+04 

-3 . 8943E+03 

1 . 0298E+02 

. N-m 

Kil 

-9 . 8025E+04 

-1 . 2847E+04 

-1 . 0298E+02 

-3 . 8943E+03' 

N-m 

Bx 

2 . 18 95E+04 

1 . 13 96E+03 

-7 . 0298E-01 

1 . 7716E+02 

N-SEC 

By 

-1 . 13 96E+03 

2 . 1895E+04 

-1 . 7716E+02 

- 7 . 0298E-01 

N-SEC 

B(J> 

-1 . 9373E-01 

6 . 8162E+01 

5 . 4617E-01 

1 . 3431E-02 

N-m-SEC 

Bl|/ 

-6 . 8162E+01 

-1 . 9373E- 01 

-1 . 3431E-02 

5 . 4617E- 01 

N-m-SEC 

Ax 

3 . 0199E+00 

-1 . 1981E- 02 

2 . 16 88E- 04 

1 . 8092E-03 

N-SEC 2 

Ay 

1 . 198 IE- 02 

3 . 0199E+00 

-1 . 8092E-03 

2 . 1688E- 04 

N-SEC 2 

A4> 

6 . 0437E- 05 

4 . 9833E-04 

3 . 5919E-05 

1 . 1575E- 06 

N-m-SEC 2 

All 

-4 . 9833E-04 

6 . 0437E-05 

-1 . 1575E-06 

3 . 5919E- 05 

N-m-SEC 2 

( CASE 

6 ) Childs 

finite length solution. 

RPM0=RPM/2 , 

l/d-i 


CYLINDRICAL SEAL 


LENGTH, DIAMETER, CLEARANCE = 1.5240E-01, 1.5240E-01, 1.9050E-04 (m) 

ROTOR, SWIRL AND DIST. SPEEDS = 3.6000E+03, 1.8000E+03, 0.0000E+00 (RPM) 

PRESSURE AT START, END AXIAL BOUNDARIES = 3.4400E+06, 0.0000E+00 (Pa) 

VISCOSITY = 1 . 2950E-03 (Pa-SEC), DENSITY = 1.0000E+03 (Kg/m3) 

FLOW = 1 . 7 71 IE - 03 (m 3 /SEC) 

TORQUE = 6 . 9241E+00 (N-m), FILM POWER LOSS = 2.6103E+03 (WATT) 


DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP . 

x (m) 

y (m) 

4> (RAD) 

(RAD) 

FORCE UNIT 

Kx 

1 . 0794E+07 

9 . 1778E+07 

-1 . 5134E + 06 

1 . 5966E + 07 

N 

Ky 

-9 . 1778E+07 

1 . 0794E+07 

-1 . 5966E + 07 

-1 . 5134E+06 

N 

Kef 

-4 . 6942E+05 

6 . 5852E + 05 

-4 . 7345E + 04 

5 . 3912E + 04 

N-m 

Kijr 

-6 . 5852E+05 

-4 . 6942E + 05 

- 5 . 3912E + 04 

-4 . 7345E + 04 

N-m 

Bx 

4 . 8718E+05 

1 . 0293E + 05 

1 . 1505E + 02 

8 . 0173E+03 

N-SEC 

By 

- 1 . 02 93E + 05 

4 . 8718E + 05 

-8 . 0173E + 03 

1 . 1505E+02 

N-SEC 

Bc}> 

- 6 . 5590E+01 

2 . 4967E+03 

2 . 8624E + 02 

5 . 1787E+01 

N-m-SEC 

Bijf 

-2 . 4967E+03 

-6 . 5590E + 01 

-5 . 1787E + 01 

2 . 8624E+02 

N-m-SEC 

Ax 

2 . 72G1E+02 

-2 . 157 9E + 0 0 

-8 . 7333E-02 

-2 . 8743E-01 

N-SEC 2 

Ay 

2 . 1579E+00 

2 . 7261E + 02 

2 . 8743E-01 

-8 . 7333E- 02 

N-SEC 2 

Acj) 

4 . 7908E-02 

1 . 6434E-01 

1 . 3703E-01 

-1 . 7189E-03 

N-m-SEC 2 

Ai|j 

- 1 . 6434E- 01 

4 . 7908E-02 

1 . 7189E-03 

1 . 3703E-01 

N-m-SEC 2 
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{ CASE 7 ) Childs finite length solution, RPM0=0, L/D= . 2 
CYLINDRICAL SEAL 

LENGTH, DIAMETER, CLEARANCE = 3.0480E-02, 1.5240E-01, 1.9050E-04 (m) 

ROTOR, SWIRL AND DIST. SPEEDS = 3.6000E+03, O.OOOOE+OO, O.OOOOE+OO (RPM) 

PRESSURE AT START, END AXIAL BOUNDARIES = 3.4400E+06, O.OOOOE+OO (Pa) 


VISCOSITY = 1 . 2950E- 03 (Pa-SEC) , DENSITY = 1.0000E+03 (Kg/m3) 


FLOW = 

3 . 9890E- 03 

(m 3 /SEC) 




TORQUE 

= 3 . 6677E+00 

(N-m) , 

FILM POWER 

LOSS = 1 . 3 827E + 03 (WATT) 


DYNAMIC COEFFICIENTS ( 

FORCE UNIT / 

DISP. UNIT ) 


DISP. 

x (m) 

y (m) 

<(> (RAD) 

f (RAD) 

FORCE UNIT 

Kx 

1 . 8583E+07 

-3 . 0274E+05 

-3 . 0756E+04 

1 . 4106E+06 

N 

Ky 

3 . 0274E+05 

1 . 8583E+07 

-1 . 4106E+06 

-3 . 0756E+04. 

N 

K0 

2 . 1014E+03 

9 . 6616E+04 

-3 . 8676E+03 

-1.6999E+01 . 

N-m 

Ki|r 

-9.6616E+04 

2 . 1014E+03 

1 . 6 999E+01 

-3 . 8676E+03 

N-m 

Bx 

2 . 1892E+04 

8 . 5184E+02 

1 . 6580E+00 

1 . 7615E+02 

N-SEC 

By 

-8 . 5184E+02 

2 . 1892E+04 

-1. 7615E+02 

1 . 658 0E+00 

N-SEC 

B0 

-1 . 0540E+00 

6 . 8066E+01 

5 . 447 9E- 01 

1 . 9870E-02 

N-m-SEC 

B\\! 

-6 . 8066E+01 

-1 . 0540E+00 

-1 . 9870E- 02 

5 . 4479E- 01 

N-m-SEC 

Ax 

3 . 0025E+00 

-8 . 3743E-02 

-5.4 944E- 04 

1 . 8925E-03 

N-SEC 2 

Ay 

8 . 3743E- 02 

3 . 0025E + 00 

-1 . 8925E-03 

-5.4 944E - 04 

N-SEC 2 

A0 

3 . 7677E-04 

4 . 5070E-04 

3 . 5698E-05 

-2 . 5857E- 06 

N-m-SEC 2 

Aijr 

-4 . 5070E- 04 

3 . 7677E-04 

2 . 5857E- 06 

3 . 5698E-05 

N-m-SEC 2 

( CASE 

8 ) Childs 

finite length solution, 

RPM0=0 , L/D=l 



CYLINDRICAL SEAL 


LENGTH, DIAMETER, CLEARANCE = 1.5240E-01, 1.5240E-01, 1.9050E-04 (m) 

ROTOR, SWIRL AND DIST. SPEEDS = 3.6000E+03, O.OOOOE+OO, 0.0000E+00 (RPM) 

PRESSURE AT START, END AXIAL BOUNDARIES - 3.4400E+06, O.OOOOE+OO (Pa) 

VISCOSITY = 1 . 2 950E- 03 (Pa-SEC), DENSITY = 1.0000E+03 (Kg/m3) 

FLOW = 1 . 7673E- 03 (m 3 /SEC) 

TORQUE = 7 . 8249E+00 (N-m) , FILM POWER LOSS = 2.9499E+03 (WATT) 


DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT ) 


DISP . 

X (m) 

y (m) 

<j> (RAD) 

l|f (RAD) 

FORCE UNIT 

Kx 

1 . 3250E+07 

7 . 5182E+07 

-1 . 8 8 82E+06 

1 . 5946E+07 

N 

Ky 

-7 . 5182E+07 

1 . 3250E+07 

-1 . 5946E+07 

-1 . 8882E+06 

N 

K<t> 

-3 . 3283E+04 

7 . 0821E+05 

-4 . 82 10E+04 

4 . 1331E+04 

N-m 

K V 

-7 . 0821E+05 

-3 . 3283E+04 

-4 . 1331E+04 

-4 . 8210E + 04 

N-m 

Bx 

4 . 8952E+05 

8 . 9262E+04 

- 1 . 4507E+02 

7 . 9030E + 03 

N-SEC 

By 

- 8 . 9262E + 04 

4 . 8952E + 05 

- 7 . 9030E+03 

-1 . 4507E+02 

N-SEC 

B$> 

2 . 5693E+02 

2 . 5639E + 03 

2 . 8825E + 02 

4 . 4699E + 01 

N-m-SEC 

B V 

-2 . 5639E+03 

2 . 56 93E + 02 

-4 . 4699E + 01 

2 . 8825E + 02 

N-m-SEC 

Ax 

2 . 72 14E+02 

-3 . 2538E + 00 

-1 . 5739E-01 

-2 . 5660E-01 

N-SEC 2 

Ay 

3 . 2538E+00 

2 . 7214E + 02 

2 . 5660E-01 

-1 . 5739E-01 

N-SEC 2 

Ad 

6 . 4178E-02 

1 . 5751E-01 

1 . 3636E-01 

-3 . 2263E-03 

N-m-SEC 2 

Ad 

- 1 . 57 5 IE- 01 

6 . 4178E- 02 

3.2263E-03 

1 . 3636E-01 

N-m-SEC 2 
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